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FOREWORD 

Man's attempts to cope with the increasing complexity of our 
present society have been more and more dependent upon electronics. 
Through sophisticated computing machines, more efficient trans- 
mission methods and automation, electronics is providing logical 
extensions of man's innate capability to process and control informa- 
tion. Implicit in these expanding relationships is the need for more 
efficient communication between man and machine. 

Computer-aided circuit design — the topic of this seminar — is 
representative of such a logical extension. The variety of newly devel- 
oped processing techniques capable of providing large-scale, complex 
electronic systems necessitates even higher levels of sophistication in 
systems organization and logical design if the opportunities for improve- 
ments in reliability, size, maintainability, and cost are to be realized. 

We are most appreciative of the favorable response on the part 
of the individual contributors to this seminar, each of whom was 
selected on the basis of important contributions to this aspect of 
computer technology. Summaries of the presentation material have 
been compiled in this document. 

We are grateful to all those present, both for indication of 
interest in the proceedings and for the evidence of active participation 
in this growing field. __ 
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Cambridge, Mass. 

James Dobbie 

General Electric Computer Division, Phoenix, Ariz. 
William Happ 
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COMPUTER METHODS OF NETWORK ANALYSIS 

By 

DR. FRANKLIN H. BRAN IN, JR. 

Dr. Branin received his B.A. degree in Chemistry from 
Wesleyan University in 1942. After several years as a radar- 
electronics officer in the Navy, he returned to graduate school 
and received his M. A. and Ph. D. degrees in Physical Chemistry 
from Princeton University in 1948 and 1950. From 1948 to 1950, 
he also worked at RCA laboratories on transistor technology and 
was recalled to active duty at the Naval Ordnance Laboratory, 
White Oak, Maryland during 1951. From 1952 to 1957, he worked 
first at Los Alamos Scientific Laboratory and then at Shell 
Development Company, Emeryville, California, on electronic 
instrumentation. 

Since 1957, he has worked at IBM, Poughkeepsie, New 
York, mainly on computer methods of network analysis. His 
contributions to this field include papers on Kron's method, 
a network approach to structural analysis, a new method of 
A-C analysis, and a network-topological description of the 
vector calculus and Maxwell's equations. 

Dr. Branin is a member of Phi Beta Kappa, Sigma Xi, 
RESA, and a senior member of IEEE. 
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This is a tutorial paper * which highlights 
the influence of computer s not only on the modus 
operand! of circuit design but also on network 
theory. The impact of computers on circuit design 
techniques is obvious enough from the proliferation 
of network analysis programs. But the influence 
of this new technology on network theory itself is 
perhaps more subtle. 

In one direction, the computer is helping to 
popularize the use of matrix notation in network 
theory because this notation is valuable in pro- 
gramming and inherently elegant for describing 
network problems. In another direction, the 
computer is hastening the recognition of the broad 
scope of network theory as exemplified by programs 
employing a network approach to mechanical and 
structural analysis. Finally, the computer is 
forcing the development of analytical and numerical 
methods which are particularly well adapted to 
both the requirements and the capabilities of auto- 
matic digital computation. 

This paper stresses the importance of keeping 
the peculiarities of the computer strongly in mind 
when choosing techniques and methods of analysis 
instead of blindly using the computer to implement 
traditional methods. Only in this way can the real 
power of the computer be properly matched to the 
problems which it is called upon to solve for the 
circuit designer. 

In keeping with this philosophy, a brief review 
is given of a matrix-topological formulation of the 
network problem which has been used by such pro- 
grams as NET-1, ECAP, and STRESS -- the latter 
having been developed for structural analysis. 

This formulation, which is applicable to networks 
of any size and complexity, is based on the pioneer- 
ing work of Gabriel Kron and on the more recent 
contributions of Roth. 

From a physical point of view, the network 
problem is concerned with predicting the behavior 
of a system of interconnected elements in terms of 
the element characteristics and the manner of inter- 
connection of these elements. Mathematically, the 
network problem involves the properties of an under- 
lying topological structure, or linear graph, and a 
superimposed algebraic structure consisting of the 
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interrelations between the network variables. The "" 
matrix approach nicely correlates these two points 
of view and precisely identifies the separate roles 
of the element characteristics and of the inter- 
connections in determining the overall behavior of 
the network. 

The bookkeeping chores associated with the 
interconnection properties of the network, and 
expressed as Kirchhoff's laws, are completely 
and automatically handled by the use of four topo- 
logical matrices: the branch-node connection 
matrix A, the node-to-datum path matrix B T (for 
the tree of the network), the circuit matrix C, and 
the cutset matrix D. These matrices are intimately 
related to one another and can be derived in se- 
quence by easily programmed algorithms from a 
table of branch node connections. 

Ohm’s law may also be expressed in matrix 
form wherin the off-diagonal matrix terms account 
for the so-called "dependent sources" that are 
frequently used to describe the behavior of active 
devices. Thus active networks are as readily handled 
as passive ones. The interrelations between the 
variables involved in expressing Ohm's law and 
Kirchhoff's laws comprise the algebraic structure 
of the network problem and are neatly summarized 
by a transformation diagram developed by Roth. 

The analysis of the electrical network problem 
by the classical mesh, node, and cutset methods 
is reviewed next. This is followed by an explana- 
tion of the mixed method of analysis which under- 
lies both Bashkow's A-matrix treatment of the 
transient problem and the so-called "state variable" 
approach. It is shown that the mixed method is 
applicable to d-c and a-c network problems as 
well as to transient problems. Thus, there are 
four general methods of analysis -- or formulation. 
Moreover, all of them can be programmed so as 
to derive the network equations automatically from 
simple input data describing the network parameters 
and configuration. 

Numerical techniques for solving both linear 
and nonlinear d-c problems are discussed briefly. 

In particular, the work of Sato and Tinney is cited 
as an example of an advance in the field of electrical 
power system analysis which electronic engineers 
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have overlooked. These authors, being forced to 
handle problems an order of magnitude larger than 
the largest usually encountered in electronics, have 
developed techniques for taking full advantage of the 
sparsity of the nodal admittance matrix. As a 
result, they can solve huge systems of equations 
(1000 or more variables) in short order and without 
exceeding the core storage capacity of a typical 
large scale computer. 

The utility of KronVs method of interconnecting 
solutions is pointed out in connection with parameter 
studies and a simplified description of the method 
is given. Although Kron's method can also be used 
to derive sensitivity formulas, an easier approach 
based on the differentiation of the nodal equations 
is presented. Finally, the solution of nonlinear 
problems by the Newton-Raphson iteration method 
and by a method developed by Katzenelson is dis- 
cussed. 

Anew approach to the solution of a-c problems 
is described next. Using the equations derived from 
the mixed method, the inverse matrix (sU - A)“* 
is computed in terms of the eigenvalues and eigen- 
vectors of the (Bashkow) A matrix, s being the 
complex frequency parameter and U a unit matrix. 
Solution of the eigenproblem for the matrix A de- 
fines a transformation which diagonalizes not only 
A but also the matrix (sU - A). As a consequence, 
the problem of computing the inverse of (sU - A) 
is reduced to the trivial task of inverting the diag- 
onal matrix (sU - A), where Ais the (diagonal) 
matrix of eigenvalues of A, followed by some 
matrix multiplications of simple form. 

This approach leads to a partial fraction expan- 
sion for any of the state variables or network func- 
tions. It permits the straightforward computation 
of the sensitivities of these variables or network 
functions with respect to frequency or any network 
parameter. Finally, the sensitivities of the net- 
work poles (which are identical with the eigen- 
values of A) can be computed. 

Since this new method of a-c analysis can be 
extended to the treatment of linear transient prob- 
lems and yields essentially the same theoretical 
information as the traditional Laplace transform 
approach, it is recommended as an alternative 
thereto, particularly for implementation on a com- 
puter. The Laplace transform technique suffers 
from two serious disadvantages when applied to 
computer- sized problems. First, the job of com- 
puting the coefficients of the polynomials in a net- 
work dunction P(s)/Q(s) is prone to serious 
inaccuracies, especially when the polynomials are 
of high degree, and may be very time-consuming. 
(For example, the "topological formula" technique 


for computing these coefficients may involve 
identifying and calculating admittance products of 
millions of trees for networks of even modest size.) 
Second, the roots of the polynomials themselves 
may be extremely sensitive to errors in the coef- 
ficients. 

By contrast, the calculation of eigenvalues and 
eigenvectors, particularly by the Q-R algorithm 
of Francis, is economical (about 10 times the cost 
of inverting a matrix) and computationally stable, 
except possibly for matrices with a very wide 
spread of eigenvalues. Moreover, the Q-R method 
readily handles the problem of multiple eigenvalues. 
This is not to say that the eigenvalue approach is 
entirely free of difficulties. But it appears to be 
far better adapted to the computer and faster than 
the Laplace transform technique. 

After a description of the analytic solution of 
the linear transient problem in terms of eigen- 
values and eigenvectors, an explanation of the 
principal obstacle encountered in the numerical 
integration of differential equations is given. It 
is pointed out that the usual integration schemes, 
such as Runge-Kutta and the predictor -corrector 
methods, are limited by the fact that the integration 
step size must be kept comparable to the smallest 
time constant (which is the reciprocal of the largest 
eigenvalue) of the network in order to avoid nu- 
merical instability. 

This limitation can be severe in networks 
having a large ratio of largest to smallest eigen- 
value. For the largest eigenvalue (smallest time 
constant) controls the permissible integration step 
size whereas the smallest eigenvalues (largest 
time constants) are principally responsible for the 
behavior of the network and determine the interval 
over which the integration must be carried out. 
Although there are some implicit integration 
formulas which are free of this restriction, they 
are of low order accuracy and require the solution 
of a system of linear equations at each time step. 
Hence, they are considered unsatisfactory for 
general use. 

A quasi-analytic approach to this problem is 
discussed in which the diagonal terms of the A 
matrix are handled analytically and the off-diagonal 
terms numerically. The stability of this method 
is governed by the largest eigenvalue of the off- 
diagonal part of the A matrix -- and this eigenvalue 
may be appreciably smaller than the largest eigen- 
value of A. Even so, by means of a similarity trans- 
formation on A, the largest eigenvalue of the off- 
diagonal part of the transformed matrix may be 
made arbitrarily small. Thus, thp integration step 
size can be increased significantly without producing 
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numerical instability. 

The effectiveness of this technique is illustrated, 
by a simple example but the problem r emains to 
find an economical method of determining and apply- 
ing the appropriate transformation(s), particularly 
in the case of nonlinear networks. 

In any event, the route that will eventually 
permit the small time-constant barrier to be breached 
appears to lie in the direction of performing mean- 
ingful operations on the A matrix itself rather than 
accepting this matrix as given and trying to devise 
clever integration formulas. For this matrix clearly 
holds the key to the entire integration problem. 

The paper concludes by pointing out some of 
the directions in which the next generation of net- 
work analysis programs should advance such as: 
incorporating improved analytical and numerical 
techniques, providing facilities for statistical 
analysis, storage and retrieval of device models, 
and the handling of much larger networks -- 


including transmission lines. It is also acknowledged 
that the advent of integrated circuits is already 
placing strong demands on existing analytical, nu- 
merical, and programming techniques and will 
inevitably force the development of more powerful 
methods for solving partial differential equations. 

Finally, the desirability of evolving from the 
present stage of network analysis to the next stage 
of optimization is recognized, but it is concluded 
that the engineer himself will continue to play an 
essential role in the design process for some time 
to come. Although the computer will certainly 
replace the routine computational tasks of the engi- 
neer, it will enhance rather than supplant his 
creative abilities. 

Reference 

1. Presented at the IEEE International Convention, 
March 20-23, 1967, in New York City. Avail- 
able from the author as IBM Technical Report 
00. 1562. 
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OLCA: AN ON-LINE CIRCUIT ANALYSIS SYSTEM 

By 

DR. HING C. SO 

Dr. So was born in Canton, China. He received the B. S. , 

M. S. , and Ph. D. degrees in Electrical Engineering from the 
University of Illinois in 1956, 1957, and I960, respectively. 

From I960 to 1965 he was on the faculty of the Department 
of Electrical Engineering at the University of Rochester. During 
the academic year 1963-64 he was on leave from the University of 
Rochester and was at the Bell Telephone Laboratories, Murray Hill, 

N. J. 


In September 1965 he left an appointment of Associate 
Professor at Rochester University to join the Bell Telephone 
Laboratories at Murray Hill, where he has been a member of the 
technical staff in the Information Processing Research Department. 
His current interest lies in computer applications and network 
theory. 
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A software system for on-line 
analysis of linear networks is described . 
Inputs to the system such as circuit 
diagrams are given via a light pen and 
a teletypewriter. Analysis results are 
available in various forms, including 
oscilloscope displays of frequency 
characteristics. The system guides the 
user through necessary steps by displayed 
comments; it also checks for user's 
errors . 

Introduction 

The application of the digital 
computer to network design is of 
widespread interest. However, current 
applications primarily use the computer 
in over-the-counter runs.l Thus input 
data are given on tapes or punched cards, 
usually according to rigidly specified 
formats. Furthermore, computation 
results from a circuit analysis program 
are available to the network designer 
only after an always-present delay, i.e., 
the turn-around time. To make the 
computer a genuinely effective tool in 
the iterative network design process, 
it is imperative that the coupling 
between the designer and the computer 
be as efficient as possible. This paper 
describes OLCA, a system of programs 
designed for on-line analysis of linear 
networks, which has existed at the Bell 
Telephone Laboratories since August 1966. 

Typically, a user of OLCA draws the 
circuit under study with a light pen on 
the cathode ray tube of a special remote 
graphic console known as GRAPHIC 1. 
Numerical data such as circuit element 
values and frequencies are supplied 
through the teletypewriter of the console. 
If desired, trial circuits representing 
deviations from the initial circuit may 
be simultaneously specified. After the 
entire problem has been defined at the 
console, it is sent to the central 
computer via data links for analysis. 

At his option the user may obtain 
displays of frequency characteristics on 
the oscilloscope of the console and/or 
tables and plots through the conventional 
output facilities of the central computer. 
Thus a user at the console can maintain 


a continuous dialogue with the main 
computer, utilizing the central computing 
power rapidly and easily to Implement his 
design. 

A full description of the experi- 
mental console, GRAPHIC 1, has been given 
elsewhere. 2 The operating environment 
for OLCA is given in Fig. 1. The central 
computer used is an IBM 7094. The 
outstanding characteristic of GRAPHIC 1 
is that it contains a small local 
computer, a PDP5, which is used to direct 
all problem-defining activities at the 
console. The central computer Is called 
into service only when circuit analysis 
is performed. This computing system Is 
in contrast to that used by the two 
on-line circuit analysis programs reported 
at MIT, 3* 4 where the Project MAC central 
computer is used to serve all console 
activities. The obvious advantage of 
not having to interrupt the central 
computer for every local task has many 
ramifications, which OLCA was designed to 
exploit fully. As a consequence, OLCA 
has the following features: 

(1) No control button board is used. 
Instead, words, called light buttons, 

are displayed on the oscilloscope, as 
shown in Fig. 2. Whenever a light button 
is pointed at by the light pen, programs 
written for the local computer enable 
the user to perform the corresponding 
function, for example, to delete unwanted 
material. Thus the user can focus his 
attention on the oscilloscope,* operation 
proceeds smoothly, and many accidental 
user errors that otherwise may occur are 
avoided . 

(2) At any stage of the problem- 
defining process, OLCA displays only 
those light buttons that the user needs. 

In addition, at any time if any action is 
expected of the user, a guiding comment 
is displayed. The user cannot proceed 
without heeding the instruction. The 
result is that the user is steered step 
by step through all operations. No 
complicated operating rules need to be 
memorized by the user. 
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Figure l.-The operating environment 
of OLCA 
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Figure 2. -A circuit being constructed. 

The words across the top of 
the oscilloscope are light 
buttons . 
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Figure 3. -A completed circuit. In this 
example the three-terminal 
device is an operational 
amplifier. 



Figure 4. -The magnitude characteristics 
of the RC -active low- pass 
filter of Figure 3 and a 
deviation circuit, as seen 
on the oscilloscope. 
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(3) OLCA checks for most accidental 
user errors locally. For example, if at 
some stage the user is expected to select 
a node, the local computer checks the 
displayed item selected to prevent the 
accidental selection of a wire or a 
component by the light pen. 

Because of these features in OLCA, 
a user with little knowledge of the 
system, and indeed with little experience 
in computing, can learn to use the system 
easily. The aim of providing efficient 
and easy interaction between the network 
designer and the central computer is 
accomplished. 

It is evident from the preceding 
discussion that OLCA consists of three 
parts : 

(a) A problem- defining program 
written for the PDP5 computer of GRAPHIC 1; 

(b) A control program written for 
the IBM 7094 for data conversion, command 
execution and result display; 

(c) A circuit analysis program5 
written for the IHyl 7094 to perform 
frequency analysis of networks. A brief 
description of these programs is given. 

A discussion on the data structure is 
also presented. 
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SUMMARY 

The design and fabrication of mask sets for 
integrated circuit manufacture has become a 
critical problem. A set of techniques utilizing 
a computer for the producing of such masks is 
described. One path is provided for the masks 
required to fabricate the devices and circuits. 
Other techniques have been developed to design 
and fabricate the wiring pattern used to inter- 
connect the circuits into a functional array. One 
such technique employes the Programmed Inter- 
connect Process (discretionary wiring) and is 
used when high levels of integration are desired 
for chips having relatively low individual circuit 
yields. 

By the use of these techniques it is possible 
to design a system which will accept circuit 
schematics as input and automatically produce 
all the masks needed to fabricate integrated 
circuit chips in order to implement the given 
circuit structure in hardware. Furthermore, it 
is possible to incorporate this automatic scheme 
into a larger system which would include such 
facilities as circuit analysis and simulation. 

INTRODUCTION 

The exploitation of integrated circuit 
technology to implement electronic hardware 
faces many challenges. One of the most serious 
obstacles is the artwork needed to generate masks 
required for the fabrication of these complex 
semiconductor structures. It is certain that with 
integrated circuit technology, we will have to 
design and fabricate an enormous number of part 
types. It is also certain that chips will be pro- 
duced with higher and higher levels of integration. 
The masks will thus be more complex and the 
likelihood of error in the design or fabrication 
of the mask sets is very great if manual tech- 
niques alone are employed. Moreover, there 
has been a growing need to decrease the turn- 
around time required from the initial design 
request for a new part until its actual fabrication. 

While it has been customary for the art- 
work to be designed by specialists who may do 
little else, there is a growing trend today in the 
need for such designs to be generated by inex- 
perienced personnel. Finally, there has been a 
growing sophistication in circuit analysis, design 
and fabrication by automatic techniques. Manual 
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mask design and fabrication represent obstacles “ 
to a unified system for the design, analysis, and 
fabrication of digital and analog componentry 
by completely automatic techniques. 

Each of these considerations argues for a 
computer-assisted process for artwork design 
and fabrication. Taken together, the need is 
overwhelming. Two sets of techniques were 
developed. One for designing and fabricating 
the artwork and the other to design and fabricate 
the interconnections between the circuits. Let 
us consider each of these techniques in turn. 

DEVICE AND CIRCUIT MASKS 

The automatic artwork generating technique 
contains three parts: A language to permit the 
designer to describe the patterns and structures 
required in the individual masks in the mask set 
to be used to fabricate devices and circuits; a 
computer procedure for translating the language 
into commands and finally, the automatic equip- 
ment which takes the commands and generates 
the masks themselves. 

The language itself is a convenient short- 
hand and symbolic way of describing the patterns 
a designer wishes to generate with a minimum 
likelihood of error. The language is easily 
learned but extremely powerful. Inexperienced 
personnel, that is, individuals who are not 
familiar either with designing masks or program- 
ming have learned the language within two days 
and have produced successful masks within 
another two days. 

The central feature of the language is a 
facility for defining individual structures and 
groups of structures on the mask and assigning 
a unique symbolic name to each group. The 
system is very analogous to that found in pro- 
gramming where a programmer defines his 
variables with symbolic mnemonic names. The 
programmer can then manipulate these names 
and leave to the computer the task of relating 
the names to the actual items of data. The 
language permits any structure or groups of 
structures to be defined by means of sub-struc- 
tures. One defines a structure which will 
generate the base of a transistor. This structure 
is then used as part of the set of structures 
defining the transistor itself. The transistor 
structures in turn can be used to define a still 
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larger structure, say, a logical decision gate. 

The transistor definition can be used as parts of 
the definitions of other structures which require 
such a type of transistor. This hierarchy of 
definitions may;bq continued indefinitely. Such 
a facility permits the gradual accumulation of a 
library of parts which can be used over and over 
again in different mask structures. Once a given 
structure has been defined it need not ever be 
redefined. Thus, in time, masks can be designed 
by appropriately choosing pieces of masks which 
have been previously designed and have been 
shown to work successfully in another and perhaps 
unrelated application. 

In general, any device which one wishes to 
construct in integrated circuit technology will 
require a number of masks . Another feature of 
the language permits any structure or group of 
structures to be represented "three -dimension - 
ally. " When a designer requests a given pattern 
to be located at a particular part of a chip layout, 
the language automatically accounts for the multi- 
plicity of masks required. The method by which 
this feature is implemented will be discussed 
below. 

The language has a replication feature which 
permits a given pattern to be repeated at a fixed 
interval around the layout. The replication can 
be either a row or column vector or a matrix 
array. Since any pattern or group of patterns is 
defined relative to an arbitrary origin, the 
replication is effected by adding multiples of the 
specified increment to the coordinate specified 
in the pattern definition. 

To initiate the process of designing the masks, 
an accurate representation of the layout is first 
drawn. With a little practice on the designer's 
part, the drawing need not be highly accurate but 
may be little more than a rough sketch. The 
designer then writes down the set of statements 
needed to produce the layout. This set may con- 
sist largely of previously prepared statements 
from other mask designs. The set of statements 
is punched on cards and fed to a computer which 
produces a digitally plotted representation of the 
masks which would be generated by the input 
statements. If any errors are found in the digital 
plot, the statements are corrected and a new plot 
of the mask is obtained. This correction proce- 
dure is repeated until the design is error -free. 

At this stage, the designer may call on the 
computer to produce a magnetic tape containing 
the commands necessary to generate the mask 
set. As mentioned before, each statement or 
statement group may refer to several masks 
simultaneously. One task of the computer is to 
remember on which mask layer the pattern is to 


be found and then sort the entire set of patterns 
mask by mask. Individual masks will then be 
produced with the appropriate patterns on them. 

The actual generation of the mask is per- 
formed by an automatic exposure device. In our 
case we have used a device called a light table. 
The apparatus is similar to an ordinary micro- 
scope. (Figure 1) A high resolution photo- 
graphic plate is placed on what corresponds to 
the stage of the microscope. A light spot is 
imaged on the plate by projecting it down the 
barrel of the microscope. The spot is turned on 
and off by an electronically controlled shutter. 
The mechanical stage is capable of motion in 
the X and Y directions and is controlled by a 
pair of digital stepping motors. The commands 
generated by the computer are decoded and fed 
to the motors and the shutter to draw the re- 
quired pattern on the photographic plate. In our 
experiment, the masks are drawn ten times size, 
the resulting plate is then used in a step-and- 
repeat camera which reduces the image a factor 
of ten and repeats the mask pattern so that a 
semiconductor wafer containing identical chips 
may be produced. 

There are clearly many advantages to this 
system. The designer can work in what is for 
him a more natural environment while the 
computer assumes the responsibility for the 
tedious bookkeeping and numerical details where 
most of the errors are made. The ability to 
gradually build a library containing previous 
designs reduces the need for re-inventing 
devices and circuits while considerably speeding 
up the design process. The replication feature 
requires the designer to furnish a minimum of 
information, further reducing the possibility of 
error. 

WIRING PATTERN GENERATION 

The method of generating the mask required 
to produce a metallized interconnection pattern 
to connect the circuits into functional arrays 
may be handled somewhat differently than in the 
case of device and circuit masks. One reason 
is the greater capability for completely auto- 
mating the design of these interconnection masks. 
Furthermore, if we wish to obtain the highest 
possible level of integration on the chip, a 
technique for generating the interconnection 
pattern in the presence of faulty circuit areas is 
required. The process that was developed to 
cope with the latter circumstance has been called 
the Programmed Interconnection Process (PIP) . 

It was recognized that the current yields of 
semiconductor devices and circuits do not per- 
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mit the interconnection of very many circuits 
directly on the chip to produce high levels of 
integration. In the PIP process, sometimes 
called discretionary wiring, a semiconductor 
wafer is fabricated with an array of circuits but 
no wiring personality. Each individual circuit is 
tested by a probe technique and the results of the 
test are fed to a computer which then determines 
which circuits are satisfactory and which are not. 
The computer then designs an interconnection 
pattern which wires together only the usable 
circuits. The output of the computer is a mag- 
netic tape which controls a light table as before. 

In this case, however, the semiconductor wafer, 
having been coated with a metallization layer and 
then a photoresist layer, is placed on the stage 
of the light table. The photoresist is exposed 
directly, hence no mask is needed to produce the 
interconnection pattern. In this process every 
wafer will have its own interconnection pattern 
determined by its particular pattern of good and 
bad circuit areas. The technique assures us that 
only functioning circuits will be interconnected to 
produce a functioning array. 

In order for the PIP process to be competi- 
tive on a cost basis, the costs unique to the PIP 
process, that is, testing the circuits, using the 
computer to generate an interconnection pattern 
and actually fabricating that pattern must be 
compatible with the other manufacturing costs. 

A major fraction of these unique costs is, of 
course, the computer time. Our PIP program 
runs on an IBM 7094 and provides interconnection 
patterns at a rate of about 0.5 sec per circuit 
interconnected. Furthermore, the program must 
be flexible enough to provide an interconnection 
pattern for every wafer, or almost every wafer, 
independent of its particular pattern of good and 
bad circuit areas. Since the PIP program must 
also be very fast, it is likely to be rather ineffi- 
cient in its use of available interconnection area. 
The requirement that almost every wafer be 
wirable further increases the amount of inter- 
connection area which is needed on the wafer. 
Hence, it is likely that a wafer layout intended for 
PIP will have a high ratio of interconnection area 
to circuit area. If we are not interested in very 
high levels of integration, or if our circuit yield 
is reasonably high, then one can assume all of the 
circuit areas in an array are good and generate 
an interconnection pattern which is the same for 
every wafer of a particular part type. One can 
then simply discard chips which are inoperative 
due to non-functioning circuits. This has been 
called the Fixed Interconnection Process (FIP). 
Here one may have the computer spend a relatively 
long time for each part type to produce a highly 
optimized wiring layout. The commands for 
generating the results of a wiring pattern can be 


run on the light table to produce a mask in a 
similar manner to the circuit and device masks 
except that the layout language need not be used 
for the interconnection mask as long as a program 
exists for generating the wiring pattern automat- 
ically. 

THE COMPLETELY AUTOMATED SYSTEM 


By combining either the PIP or FIP wiring 
approach with the library facility of the layout 
language, a completely automatic system for the 
generation of artworkfor all the masks required 
is obtained. The input to this system would only 
be the circuit schematics to be implemented. To 
visualize the process we refer to Figure 2. A 
program can be written to assign circuit types 
VERTICAL VERTICAL 

CHANNEL CHANNEL 



Figure 2. - Schematic view of an integrated 
circuit wafer. Circuit areas are 
fabricated in a rectangular array while 
the interconnections are in vertical and 
horizontal channels. 


to the circuit areas shown in the figure. Since we 
are given the circuit structure to be implemented, 
this assignment can be performed to make most 
efficient use of the layout. We can, for example, 
minimize interconnection distance between circuits 
or improve the wirability of the layout by avoiding 
wiring congestion. Knowing what circuit type is 
to be located in each of the circuit areas, we can 
use the layout language to call out a previously 
designed group of patterns producing that circuit 
type and locate the pattern, in the given circuit 
area. This is enough information to allow us to 
construct the masks which will fabricate that 
particular wafer. The interconnection pattern 
may be generated using either the FIP or PIP 
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program discussed above. The entire process 
from circuit schematic input to finished masks 
can take place entirely within the computer 
without any manual intervention being required. 
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SUMMARY 

The simulation of dynamic systems 
and continuous filter networks and the 
filtering or processing of data signals 
by means of digital computers require 
both the design and utilization of digital 
filters.* In simulation work the digital 
filters are designed to represent or 
approximate continuous systems such that 
the behavior of the continuous systems 
can be evaluated, either numerically or 
subjectively, for a wide range of input 
signals. This is in contrast to the 
more usual network transient response 
calculations based on very simple inputs 
such as steps, ramps, impulses, and 
single sinusoids. 

The effective use of the discrete 
models depends Q n a detailed knowledge 
of the nature of the discretizing 
approximation process used, of the 
sampling-frequency vs. simulation- 
accuracy tradeoff, and of the essential 
linearity or nonlinearity of the con- 
tinuous systems being modeled. This ' 
paper is, therefore, concerned with a 
brief discussion of two of the more 
commonly used digital filter design 
methods. This is followed by a general 
description of a digital filter design 
program which is capable of generating 
discrete representations for linear 
continuous filters including a wide 
variety of the classic filter forms. 

The paper concludes with comments on the 
application of the program. 

Introduction 

Digitalizations may be sought for 
both linear and nonlinear systems. For 
nonlinear systems the choice of 
sampling frequency may have to -be made 
quite large to yield sufficient accuracy 
in the simulation. Here the amount of 


*The term digital filter refers to the 
computational process or algorithm by 
which a sampled signal or sequence of 
numbers acting as an input is transformed 
into a second sequence of numbers termed 
the output signal. 


N 67-22625 

calculation required may easily approach 
that involved in directly integrating the 
differential equations by standard 
numerical techniques. On the other hand, 
the combination of the superposition 
property of linear systems and the quasi- 
band limited nature of realistic input 
signals makes possible the efficient 
simulation of systems at sampling rates 
approaching the Nyquist sampling limit. 

Linear continuous systems which are 
describable by linear differential 
equations with constant coefficients 
become representable by linear difference 
equations when in the discrete form. 

Thus no inherent additional complexity is 
added by going from the continuous domain 
to the discrete domain of representation 
[l]. Numerous techniques of digitali- 
zation have been described in the 
literature [2]. Of all these techniques, 
the two most widely used in the 
simulation of continuous filters and 
networks are the standard z transform 
and the bilinear z transform. These two 
methods differ in the nature of the 
approximation made in going from the 
continuous domain to the discrete domain. 

Design Methods 

Viewed from the frequency domain, 
the approximation made in using the 
standard z transform is easily seen from 
the defining equations 


H*(s) = ^H(s+jna> s ) + | h(0 + ) (1) 

n=-a> 

and 

H^z" 1 ) = T^h(nT)z' n (2) 

n=0 
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Thus the impulse response of the obtained 
digital filter, H*(z _1 )> is identical to 
the samples, h(nT), of the impulse 
response of the continuous filter H(s). 
However, the frequency response 
characteristic, H*(s), of the digital 
filter in general differs from that of 
the continuous filter, H(s), by the 
addition of the terms H(s+jno) s ), n/0, 
this effect, being termed folding or 
aliCsChg. Unfortunately, when H(s) is a 
rational function of s, it is not band- 
limited and therefore H(s) / H*(s) in 

/ % 

the baseband [ - < to < 1, This is 

especially true for many filter designs 
such as elliptic and Chebyshev Type II 
filters for which the high frequency 
behavior is either a constant or at most 
K/s. Application of the standard z 
transform technique consists simply of 
first obtaining a partial fraction 
expansion of H(s) in its poles and then 
z transforming each of the individual 
terms using transform pairs derived from 

1 T 

s+a 1 _ e -aT z -'l (3) 

-1 -«T 

where z - e 


H*(z 1 ) = H(s) 


2 a- z -]) 
T d-- 1 ) 


( 6 ) 


The price paid for this simplicity in 
digitalization is the nonlinear warping 
imparted to the frequency scale as can be 
seen by setting s^ = Jon and s w jco in 
(4) to yield 


<oT 

~T 



(7) 


Because of this frequency warping 
aspect (especially for frequencies greater 
than 2/T) , the bilinear z-form is most 
useful in obtaining digital filter 
approximations for continuous filters 
whose magnitude characteristics are 
essentially piecewise constant in 
successive pass and stop bands. This type 
of frequency behavior is typical of many 
lowpass, bandpass, and bandstop filters. 
Compensation can be made for the effect of 
warping by prewarping the continuous 
filter design in the opposite way so that 
upon applying the bilinear z form, the 
critical frequencies will be shifted back 
to the desired values [3]. 


A digitalizing technique which does 
not introduce error due to aliasing is 
the bilinear z transformation defined by 


tanh 




w 


which becomes upon substituting 

_-i -i T 


s 


_ 2 (1-z ) 

" T a-- 1 ) 


(5) 


Thus, In terms of the z plane, this 
algebraic transformation has the property 
of uniquely mapping the left-half of the 
s-plane into the exterior of the unit 
circle in the z~* plane. Folding errors 
are eliminated since no folding occurs. 


The bilinear z-form is applied 
simply by making the substitution 
indicated by (5) in the transfer 
characteristic H(s). Hence, 


A Filter Design Program 

The versatility of these two design 
procedures has been utilized effectively 
by Golden and Kaiser in the implementa- 
tion of a general digital filter design 
program. The program has been further 
refined by Golden to couple its results 
effectively to a simulation program for 
complex systems. The basic scheme of the 
design program is shown in Fig. 1. 



FIG.1 
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The program is divided into 
essentially two parts: the first part 

generates the linear continuous filter 
design required to meet the filter 
specifications, while the second part 
comprises the digitalization procedure 
which yields the desired coefficients 
for the resulting digital filter design. 
For input the user may specify either 
l) the actual continuous filter transfer 
characteristic whose digitalization is 
desired or 2) simply the general type 
of filter characteristic desired and the 
specifications associated with that 
filter type. The program can generate 
the basic lowpass filter types - 
Butterworth, Bessel, transitional 
(Butterworth-Thomson) , Chebyshev Type I 
and Type II, and the elliptic or Cauer 
parameter types - from specifications 
such as in-band ripple, out-of-band loss, 
transition ratio, filter type, band-edge 
frequencies, and order, if applicable. 

The program then has facility for trans- 
forming any of these lowpass designs to 
bandpass, bands top, or highpass filter 
designs by means of the standard 
frequency band transformation. The basic 
lowpass designs and band transformations 
are discussed in detail in many network 
synthesis references, for example 
Weinberg [4 ] . Provision is also included 
to prewarp automatically the original 
filter band-edge frequencies and thus 
give a prewarped continuous filter design 
if the bilinear z transform is to be used. 

The digitalization section of the 
program gives the user the choice between, 
the standard z transformation and the 
bilinear z transform with or without 
warping correction. The results of the 
digital filter design are printed out as 
the coefficients of the digital filter 
in the parallel canonical form 


H*( z” 1 ) = A oq + 


+ A 


if i + B - 


01 


li* 


This program has been a valuable aid 
in determining the digital filter co- 
efficients required for the simulation of 
complex systems via the block diagram 
compiler BLODIB [6J. When special purpose 
digital filters are to be physically 
Implemented by hardware, the program 
described can be used to obtain a first 
set of values for the filter coefficients. 
From this set of values the designer can 
vary the bit size representation of the 
coefficients until he has obtained a 
satisfactory compromise between filter 
complexity and performance characteristics. 
The requirements on filter coefficient 
accuracy have been Investigated at some 
length by Kasier [5], Other investigators 
have started to look into the problem 
associated with roundoff and truncation 
in the carrying out of the arithmetic 
operations of the digital filter. The 
overall picture of error analysis is by 
no means complete as it depends also on 
the nature of the input signal and the 
performance criterion as well as the 
realization scheme. 

In this paper the use of digital 
filters for data processing and simulation 
of complex systems has been motivated. 

Two digitalization procedures were 
briefly discussed. A digital filter 
design program was then described which 
incorporated the digitalization 
procedures. The paper concluded with 
brief comments on the range of 
application of the program. 
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The general features of the NET- 2 program 
are presented. This program, under development on 
the MANIAC computer, will perform nonlinear DC and 
transient analysis and linearized AC analysis of 
arbitrary circuit configurations. 

The language for describing the circuit sche- 
matic is presented in detail. Features for Inclu- 
sion of models of various linear and nonlinear de- 
vices are discussed. The program is written so as 
to permit device models to be added without rewrit- 
ing NET- 2 itself. Models and their parameter 
values are stored in a model library and called 
into the calculation automatically. Input facili* 
ties are included for superseding library values 
of device parameters. 


resistors, capacitors, inductors, mutual inductive 
couplings, switches, voltage and current sources, 
and a library of device models. The circuit ele- 
ments may be interconnected without restriction, 
provided a physically impossible situation does 
not result (e.g., a zero impedance mesh of voltage 
sources). 


CIRCUIT DESCRIPTION 

The circuit to be analyzed is described by a 
sequence of entries, each of which is a description 
of a single circuit element. These entries may be 
written following a process of assigning, on the 
schematic diagram, identification symbols to the 
elements and names to the connection points. 


A major feature of NET-2 is its variational 
capability. Three types of variational procedures 
are available: the STATE calculation which varies 

circuit parameters in discrete steps, the Monte 
Carlo calculation which determines the statistical 
distribution in response variables for a given cir- 
cuit parameter population distribution, and the 
Extrema Search which extremizes the value of a 
response variable within an allowable circuit 
parameter space. 

INTRODUCTION 

The NET- 2 Circuit Analysis Program is cur- 
rently under development on the MANIAC computer at 
the Los Alamos Scientific Laboratory. It is a 
logical extension of the highly successful NET-1 
program which was completed in 1962. NET-1 per- 
forms nonlinear analysis of circuits whose topol- 
ogy is fixed using nominal values of the circuit 
elements; NET-2 extends the analysis to variational 
studies in the circuit element values as well as 
permitting changes in topology. 

NET- 2 observes the same philosophy of ease of 
usage as NET-1: the user describes his circuit in 

a simple engineering- oriented language, and he is 
not required to have any knowledge of mathematical 
techniques, programming languages, or computational 
methods in order to use the program successfully. 

Many of the restrictions in NET-1 have been 
removed to permit a more natural description of the 
circuit. A circuit described to NET-2 may exist in 
several electrically isolated parts. Voltage and 
current sources which are floating in the circuit 
are handled in a straightforward manner. Further- 
more, circuits are not required to contain either 
sources or impedances, thus permitting the handling 
of switching and path problems. 


Circuit elements are identified by an ID con- 
sisting of a letter prefix, designating the kind of 
element, and a non-zero integer. The integers need 
not be sequential and the same integer may be used 
in ID for different kinds of elements. The ID may 
not contain more than six characters. 

Voltage and current sources are assigned ID in 
the form Vl6 and 14. Linear elements such as re- 
sistors, capacitors, inductors, mutual inductive 
couplings, and switches are assigned ID in the 
form R34, Cl6, L34, K5, and S2. Modelled devices, 
such as transistors and magnetic cores, are assigned 
ID in the form X12, X7, etc., where X is the appro- 
priate letter prefix for the particular device. 

All connection points are assigned arbitrary 
alphanumeric names consisting of a maximum of six 
characters each. The ground node is assigned the 
integer 0 as its name. It is not necessary that 
there be a circuit path between various isolated 
parts of a circuit. 

If a source element is connected between a 
node and ground, it is convenient to give the node 
the same name as the source, such as V12, l6, etc. 
Such a designation allows a simplification in the 
description of the source element itself. 

Parallel circuit segments exist whenever a 
circuit contains several identical segments, and 
these segments have all connection points into the 
rest of the circuit in common. NET- 2 circuit ele- 
ment description formats have special provisions 
for describing such parallel segments. In describ- 
ing the parallel segments, node names and IDs are 
assigned within only one of the segments. The 
element descriptions for that segment are then 
tagged with a number in parentheses to indicate 
how many parallel segments are represented. 


*Thls work performed under the auspices of the The circuit element formats for resistors, 

U.S. Atomic Energy Commission. 
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capacitors, and. inductors are: 

Rn (p) f g xxx 

Cn (p) f g xxx 

Ln (p) f g xxx 

where Rn, Cn, 'in j# element ID 

p » numbeir of parallel units (optional) 
f and g » connection points 
xxx = value of element in kilohms, pico- 
farads, or microhenries. 

The circuit element format for mutual induc- 
tive coupling is: 

Kn i 3 xxx 

where Kn — coefficient of coupling ID 

i and j = IDs of the inductors involved 
xxx = value of the coefficient of coupling. 

The switch is an ideal, zero impedance, single 
pole single throw device. It may he initially 
open or closed during a steady state analysis as 
well as operated at prescribed times during the 
transient analysis. Ihe format is: 

Sn f g Mode 
Mode xxx 
Mode yyy 
Mode zzz 

where Sn = switch ID 

f and g = connection points 
Mode e status of the switch, specified by 
the words "Open" and "Close." The 
first mode designation, following the 
g connection point, is the setting of 
the switch during the steady state 
analysis and at the initial time of 
the transient analysis. Subsequent 
mode designations indicate opening and 
closure of the switch at the respective 
times xxx, yyy, ;-;zz, etc., during the 
transient analysis. 

Source elements may be constant in value or 
time varying. The general format for voltage 
sources is: 

Vn f g Value 

Current sources may also include parallel segment 
information: 

In (p) f g Value 

where Vn and In = element ID 

p = number of parallel units (optional) 
f and g = connection points 
Value = information regarding waveshape and 
associated amplitude and time param- 
eters. 

The value of the source is always referred to con- 
nection point g. If the name of connection point 
f is the same as the source ID, and connection 
point g is node 0, it is permissible to omit the 


specification of f and g in the source element 
description. 

The Value portion of the source element de- 
scription specifies the time behavior of the 
source. Sources may be constant in value, or be 
time dependent. Time dependent sources include 
repetitive trapezoidal pulses, nonrepetitive 
tabulated waveforms, decaying exponentials, and 
sine waves. In addition, a small signal sine 
rave is available for excitation of the AC steady 
state analysis. Both amplitude and phase of this 
AC sine rave may be specified. 

NET-2 has the capability of utilizing modelled 
devices which may be defined at any time and 
placed in a model library. This process is a 
relatively simple one and does not necessitate 
the rewriting of the NET-2 program. The models 
may represent multi-terminal linear or nonlinear 
devices, such as voltage controlled current 
sources, silicon controlled switches, magnetic 
circuits, photosensitive devices, and entire 
functional circuit blocks. 

Each modelled device is represented in the 
analysis by an equivalent circuit, a set of con- 
trolling equations for the equivalent circuit, 
and a set of device parameters for the coeffi- 
cients of the controlling equations. The general 
format for device entries consists of the device 
ID, the number of parallel units, the connection 
points (as many as needed), a type name, and op- 
tional mode information. The ID prefix specifies 
the kind of device (e.g., klystron, nonlinear in- 
ductance, etc.) from which NET supplies an equiva- 
lent circuit and controlling equations. The type 
name enables NET to obtain particular parameter 
values for the equations so that different devices 
of the same family (e.g., 2N705, 2N706A, etc., in 
the case of transistors) may be represented in 
the calculation. The mode information specifies 
constraints which may be imposed on the device 
during analysis, e.g., the specification of the 
magnetic state of a magnetic core. 

The device parameter values may be changed as 
part of the input by means of the PARAMETER entry. 
In this entry specific devices are indicated, fol- 
lowed by the parameter name and its new value, for 
example : 

PARAMETER 


T35 

Bn 25 
Ics ,001 

D3 

Temp 35 


The circuit description provides for a nomi- 
nal nonlinear DC steady state solution and a 
transient solution using the nominal DC solution 
for initial values. If the circuit contains an 
AC sine source and a specified frequency value for 
this source, a linearized AC steady state solution 
will also be produced. This linearized AC solu- 
tion is obtained by first solving the DC nonlinear 
problem, evaluating the impedances at the DC 
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operating point, and then solving the small signal 
AC steady state problem using these Impedance 
values . 


VARIATIONAL CAPABILITIES 

It is possible to generate additional AC and 
DC steady state solutions by varying source values, 
circuit element values, device parameter values, 
and circuit topology. There are several ways of 
accomplishing the variation, depending upon the 
nature of the study undertaken. 

The STATE entry provides one method of modify- 
ing circuit values and topology. In this entry 
specific values for circuit elements and parameters 
are indicated. These new values supersede the 
nominal values. Any values not specifically en- 
tered in a STATE entry are taken from the nominal 
circuit description. 

It is also possible to vary circuit parame- 
ters and frequency in a set of discrete linear or 
logarithmic steps in a STATE entry so as to per- 
form a parametric study. It is possible to refer 
to the entire ensemble of any element kind by 
using only the letter prefix for that element. 
Changes in device parameters may be made with the 
same general format used in the PARAMETER entry. 

The order of listing all items is important since 
the last value listed is the one used. 

Device modes may be specified, but such infor- 
mation must be listed separately from any parame- 
ter changes. The word NONE may be used to remove 
a previously specified mode. 

An example of a STATE entry is: 

STATE 5 


R7 

1.5 


T 

D2 


HG 

Temp 

35 

T 

Bn 

68 

T2 




Bn 


R 

.1 w 

.5 

R5 

2.7 


T 

OFF 


T2 

NONE 


S12 

1 CLOSE 



This list is interpreted in the following manner: 
Temp = 35 for diode D2, all transistors and 
all Hall generators 
Bn = 68 for all transistors except T2 
Bn - 4-5 for T2 

The value of all resistors except R5 varies in 
4 linear steps from .1 to .5, giving a 
total of 5 solutions. All resistors 
which are varying will have the same 
value during any given solution. Note 
that R7 is not held at 1,5 because the 
variational entry for all resistors 
supersedes the R7 entry. 


R5 = 2.7 

All transistors except T2 are held off. Note 
that the NONE tag cannot be combined 
with the substitution for Bn of T2. 

S12 is closed. 

It is possible to generate two dimensional 
studies (or higher dimension studies) as in the 
example: 

STATE 3 

0V 12 -30 (if) -60 

T4-5 

Bn 75 

Ics .001 (5) .002 

In this case the phase of source V12 and the Ics 
parameter of T45 will be varied over a two dimen- 
sional grid of size 5x6, giving a total of 30 
solutions. 

If a logarithmic variation is desired the 
number of steps is suffixed with an asterisk: 

STATE 57 

FREQ .001 (25*) .1 

In this example the frequency of the AC sine 
sources is varied logarithmically from .001 to .1 
in 25 steps. 

Another variational calculation available in 
NET- 2 is the Monte Carlo calculation. Here NET 
randomly picks circuit parameter values from known 
population distributions, constructs and analyzes 
the resultant circuit, and then tabulates the pop- 
ulation distribution for the circuit responses of 
interest. The resultant circuits will then have 
response distributions which are statistically 
identical to the actual circuit in mass production. 

In order to accomplish the Monte Carlo calcu- 
lation, NET must know the population distribution 
shapes for the circuit parameters, the numerical 
values at the end points of the distributions, and 
the number of circuits which should be analyzed 
to give the representative statistical behavior. 

Distributions available include flat, trian- 
gular, Gaussian, and tabulated distributions. 
Distributions may be specified in either linear or 
logarithmic space. Limit values may be specified 
a6 either actual values or relative (percentage) 
values. Population distributions for device 
parameters may be stored as part of the model 
library. 

A third variational capability in NET- 2 is 
the Extrema Search. In this type of calculation 
NET attempts to extremize the value of a particular 
response variable by varying circuit parameters 
vithin an allowable parameter space. Since the 
extrema point for a given response variable may 
not be at the limit values of the circuit parame- 
ters, this method is more powerful- than the so- 
called worst case analysis. It is also valuable 
in the optimization of circuit response. Both DC 
and AC responses may be treated with this calcula- 
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tion. 


INPUT- OUTPUT 


TRANSIENT ANALYSIS 

Transient analysis always proceeds from the 
initial values calculated by the nominal DC analy- 
sis. The transient analysis can be terminated 
automatically when a termination condition has 
been satisfied. The termination condition may be 
composed of simple Boolean functions of a variety 
of criteria, including the switching status of de- 
vices, values of response variables, and elapsed 
real time. 


Input-output capabilities depend upon hard- 
ware availability in specific computers. The 
pilot version on the MANIAC computer has input via 
cards, paper tape, and on-line keyboard. Output 
may be in printed or graphical form. Transient 
solutions may be displayed as computation is in 
progress. Graphical displays include linear, 
semilog, log- log, and histogram plots for all re- 
sponse variables, node-to-node impedances, sensi- 
tivity coefficients, and functions of these quanti- 
ties. Variational steady state analyses may be 
displayed in parametric form. 
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SUMMARY 

A new method of sensitivity calculations 
for single parameter and multiple parameter 
sensitivities has been developed. The method 
has a direct interpretation in terms of network 
configuration. Furthermore, it may be used on 


E e (k)* = E b (k) + E s (k) 
I e ( fc) - I s (k) + I b (k) 
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k = 1 , . . . , n 

I (k). 


(la) 

(lb) 


existing digital computer network analysis pro- 
grams without modification. The ease with which 
multiple parameter sensitivity calculations can be 
made has been used to analyze continuously equiva- 
lent networks and the mathematical model of the 
human cochlea. 

INTRODUCTION 

In performing a sensitivity analysis, it is 
often desired to calculate the sensitivity func- 
tion of a large number of response variables as 
a function of a single parameter. Conversely, 
the sensitivity function of a single response 
variable to changes in a large number of para- 
meters may be desired. Efficient solutions to 
both problems are presented in this note. The 
solution to the first problem is quite general in 
applicability. The solution to the second problem 
is limited to AC calculations for reciprocal net- 
works unless one wishes to do convolutions. 



One distinct advantage of the method is 
that currently available digital network analysis 
programs, such as ECAP and NET-1, can be used 
without modification for sensitivity calculations. 

Although examples in this paper are limited 
to electrical networks, any analagous systems 
may be analyzed by the same method. The sensi- 
tivity functions in this paper are obtained with 
respect to the fundamental parameters of the 
network, resistance, capacitance, inductance, 
and current gains. 

DEVELOPMENT OF SINGLE PARAMETER SENSITIVITY [1] 

The connected electrical network to be con- 
sidered consists of N branches, each contain- 
ing a linear, time invariant resistor, inductor, 
or capacitor. The extension to nonlinear/ time 
varying networks is trivial. However, existing 
programs are not usable without modification. 

Each branch may also contain an independent 
current and/or voltage source, as well as a 
dependent current source. The configuration of 
such a standard branch is shown in Figure 1> 

The branch relations for the kth standard 
branch are 


Figure 1. -Standard branch 

The element equations for the standard branch 
will be in one of the following forms: 

E (k) = R(k) I (k) (2a) 

« (k) 

E e (k) - L(k) -jt— k . ! n <»> 

d E <k) 

I e (k) = C(k) — jj (2c) 

In addition, the controlled current source 
must be considered. It is not an element in the 
same sense as R(k) , C(k) or L(k). The equation 
for the controlled current is: 

V k > ■ Bmk < 2d > 

In order to find the voltage or current 
sensitivity of the kth branch with respect to an 
element variation in the jth branch we take par- 
tial derivatives with respect to A(j), the element 
value of the jth branch. 


dE e (k) _3E b (k) dE s (k) 
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Rewriting the above equations in more com- 
pact form, with derivatives of the independent 
sources set to zero, we have 

E'(k) = E'(k) (3a) 

e b 

i;oo i^(.k) - r<k> (3b) 

Note that I^(k) will have the form 


i;<k> 

- 

V * A(J > 

(4a) 

VU) 

= A(j)i;(m) + I e (m) 

Vi • 

(4b) 


Note that the controlled current source is con- 
sidered part of the branch relations and is not 
considered an element. 


Differentiating the element equations, we 
obtain 


original network in three ways: 

(1) . All independent sources of the original net- 

work have been set to zero. 

(2) . The jth branch now includes a new driving 

source which may be represented in one of 
the following ways, depending on the element 
type (A(j) i 4 P m j): 

a. Resistor: A voltage source of value 

I (j), placed in series with the element 

R?j)*. 

b. ^nducttjr: A voltage source of value 

— , placed in series with the ele- 

ment L(j). 

c. Capacitor: A current source of value 

d£ (j) 

— ^ , placed in parallel with the 

element C(j). 


e;(r) 

e;<j> 

E;(k) 

i;o) 

i;o> 


R(k) 

R(j) 

L(k) 


i;o) + 

d i;(k) 

dt* 


I (J) 


MJ> 

0(k) 

C(j) 


d i;(j) d I e (j) 
dt + dt 

d E^(k) 
dt 

d b;.(j) d E e o) 

dt + dt 


A( j) * R(k) 
A(j) = R( j) 
A(j) t L(k) 

A( j) = L( j) 
A(j) t C(k) 
A(j) = C(j) 


(3). If the f3 . «= A( j) , a current source of 

value I is placed in parallel with the 

existing controlled current source. 

To apply the above results to a sensitivity 
analysis, it is convenient to express the voltage 
driver in the jth branch in terms of its Norton 
equivalent. Thus, the procedure for calculating 
network sensitivities to variation in a single 
component reduces to the following: 

1. Construct an "auxiliary’' network which 
is a duplicate of the original network, 
with all independent sources set to 
zero (deleted).’ 



(5) 

We may describe the network 

in matrix form 

with the following equations 


[B)E b ] = 0 

(6a) 

[S]I b ] = 0 

(6b) 

where [B] and [S] are the circuit 

and cut-set 


matrices for the network; E, and are vectors 

of branch voltage and branch currents respective- 
ly. Differentiating (6) with respect to A(j) 
we obtain the network sensitivity equations 


[B]E£] = 0 

(7a) 

[S]I£] = 0 

(7b) 


Drive this "auxiliary" network by means 
of a dependent current source placed 
across the jth branch. The source is 
proportional to the current flowing in 
the jth branch of the original network 
if A( j) £ The constant of propor- 

tionality is - /*;. • for resistors and in- 


ductors, and ~y for capacitors, 

where A(j) is the value of the element 
in the jth branch. For (3 . = A(j), 
the current source is I e (m^ 


3. The voltages and currents in the "auxili- 
ary" network are respectively the voltage 
and current sensitivities of the network 
to a variation of the parameter A(j). 


Substituting (3a) and (3b) into (7a) and (7b), 
we have 


DEVELOPMENT OF MULTIPLE PARAMETER 
SENSITIVITY CALCULATION [2] 


[b]e;i = o 

(8a) 

tsu;i + [sir] - o 

(8b) 


Considering the element sensitivity relations 
(5) and the topological equations (8a) and (8b), 
we see that the equations for the sensitivity 
calculation differ from those describing the 


By definition, a reciprocal network has the 
following property: If an excitation is applied 

at one pair of terminals in the network, and a 
response is measured at some second pair of termi- 
nals in the network, interchanging the points of 
excitation and response does not change the ratio 
between these quantities [3], Now consider cal- 
culation of the sensitivity of the voltage across 
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branch B with respect to parameter A in branch A. 
By the reciprocity theorem, the voltage that 
appears across branch B when branch A is excited 
by the current I^/A will be the same as the 
voltage that appears across branch A when branch 
B is excited by the current I. /A (see Figure 2). 
Therefore, in the auxiliary network used to cal- 
culate the sensitivity function, placing the 
current source of value I /A in parallel with 
element B rather than in parallel with A, results 
in the branch voltage across element A being the 
voltage sensitivity of the voltage across B with 
respect to element A. 



Fig. 2. Property of Reciprocity for Element A in 
Auxiliary Network 


Consider now calculating the A-C sensitivity 
of branch B with respect to all other parameters, 

A, B, C, D,... in the network. Suppose in 

the auxiliary network, a current source of unity 
value and zero phase (abbreviated 1/0) is placed 
in parallel with element B (see Figure 3). Ob- 
viously, the voltage E* does not represent the 
previous voltage sensitivity B Eg/BA any more, 
nor is E* equal to 3Eg/BC, etc. However, the 
voltages E*, E*, etc. are related to the sensi- 
tivity of Eg with respect to A, C, etc. re- 
spectively. If the current had been I ./A instead 
of 1, then E* would have been BEg/dA. Thus, 
if E* is multiplied by I /A, the resulting 
quantity (E*) (1^/A) is equal to 3Eg/3A. 


Similarly, if E* is multiplied by Ip/C, the 
result is BEg/BC (se Figure 3). 



Fig. 3. Auxiliary Network with Unity Value 
Current Source 


Therefore, to obtain the A-C sensitivity of one 
response voltage with respect to all elements in 
a reciprocal network, one has the following rules: 

1. Place a current source of value unity with 
zero phase in the auxiliary network in paral- 
lel with the element of which the sensitivity 
of the response function is required. 

2. The voltage across any other element in the 
auxiliary network multiplied by the current 
through the same element in the original net- 
work and divided by the element value itself 
will give the voltage sensitivity. (-1 for 
capacitors) 

3. Repeating (2) with all elements the voltage 
sensitivity of a single response variable to 
all elements can be obtained. 

EFFICIENT PROGRAMMING OF SENSITIVITY 
CALCULATIONS 

In the case of single parameter sensitivity 
calculations, the method developed above consists 
of constructing an auxiliary network and driving 
this network with a controlled current source pro- 
portional to the current in the element whose 
sensitivity is desired. This interpretation pro- 
vides the means by which one can use existing pro- 
grams, such as ECAP or NET 1 without modification. 
However, it does increase the size of the network 
which must be calculated. Another way to do the 
problem is to form the sensitivity equations when 
the network equations are formed by adding dummy 
driving current sources to the original network 
across the elements whose sensitivity is required. 
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Then, this set of equations is solved with the 
dummy sources set equal to zero. Each current 
through the elements whose sensitivity is re- 
quired is stored. Then, the network is solved 
recursively with the normal driving sources set 
equal to zero but the desired sensitivity source 
not equal to zero. Thus, no network larger than 
the original network has to be solved. This 
allows one to use the full capability of the net- 
work program to solve a sensitivity problem. 

This is the method used in RAPID [4], 

APPLICATIONS 

Much literature exists on the applications 
of sensitivity calculations. Therefore, two 
applications of the sensitivity analysis method 
as presented above, will be given. The applica- 
tions are a result of the ease that one can cal- 
culate multiple parameter sensitivities. The two 
examples are, (1) the continuously equivalent 
network and its sensitivities [2] and (2) evalu- 
ation of mathematical models [5], 

The continuously equivalent network is a 
network in which the elements are varied as a 
continuous function of some dummy variable, but 
whose specified properties remain constant. By 
defining a performance index which is a function 
of the sensitivity function, the "best" continu- 
ously equivalent network can be selected out of 
all possible ones. Previous work examined this 
problem at one frequency [6], With the new tech- 
nique, these networks were examined over a range 
of frequencies. The experimental results of this 
study allows one to make three hypotheses about 
continuously equivalent networks: 

1. The continuously equivalent network result- 
ing from a minimization of the sum of the 
magnitudes squared of the sensitivity func- 
tions at a given frequency is the network 
with minimum sum of the magnitudes squared 
of the sensitivity functions at all fre- 
quencies. 

2. The sum of the magnitudes squared of the 

sensitivity functions decreases as the num- 
ber of elements increases in continuously 
equivalent networks. 

3. The sum of the sensitivity functions is in- 
variant with respect to the various equiva- 
lent networks. 

The most important hypothesis is the first. 
If true for all continuously equivalent networks, 
then a great simplification of the computation 
problem results. Even if a limited class of con- 
tinuously equivalent networks exhibit this prop- 
erty, then the result is worthwhile. These 
experimental results do not prove that all con- 
tinuously equivalent networks have these three 
properties, but they do provide a challenge for 
future work. At least they indicate a class of 
continuously equivalent networks which have this 
property. 


In recent years several investigators have 
constructed models to represent the response of 
the human cochlea [7] [8], These models are par- 
tial differential equations which, in one case at 
least, have been reduced to an electrical analog. 
As an example of what could be done with sensiti- 
vity calculations, the electrical analog system 
was solved using ECAP. Then, using the method of 
multiple sensitivity analysis, the response of 
the network was calculated and compared with ex- 
perimental measurements. Then, a linear program- 
ming optimization routine was used to adjust the 
parameters in order that the error between the 
model and the experimental measurements was mini- 
mized. The result was that the experimental 
measurements and model agree quite closely in 
magnitude, but not in phase. It is conjectured 
then, the errors in phase are not due to model 
errors but are due to measurement errors [6j. 
These results will be presented in detail else- 
where. 
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SUMMARY 

CIRCUS is a digital computer program which 
simulates the time domain response of an elec- 
tronic circuit to. an arbitrary forcing function. 
Although motivation for the development of CIRCUS 
originated from efforts to evaluate transient 
nuclear radiation effects on electronics, the 
program Is applicable to the much wider class of 
everyday electronic design problems, and any 
physical system having an electrical analogue. 

Like other time domain circuit analysis 
programs such as NET-1, PREDICT, and SCEPTRE, 
CIRCUS requires a topological description of the 
circuit which specifies element values and 
forcing functions, and from this information 
constructs the appropriate time-domain equations, 
finds the dc initial conditions, and then evalu- 
ates the transient solution by numerical 
integration of the differential equations. Sig- 
nificant features of the program are its speed, 
flexibility, easy input and output programming, 
and present availability on four different com- 
puter systems. 

Introduction 

In the past few years, the spurious effects 
generated by transient nuclear radiation on 
electronic systems have come under intensive 
study by designers of military electronic sys- 
tems. Although the radiation pulse incident 
upon a given electronic system is transient, the 
effects it causes in the electronics may be both 
transient and permanent, and these effects may 
lead to either transient or permanent system 
failure. To analyze such effects, it was evident 
that some means of establishing the vulnerability 
of electronics to radiation without extensive and 
costly system testing was required. The attempt 
to analyze complex circuitry by hand presented 
considerable difficulties. Eventually, digital 
computer aids were selected as the most practical 
way to solve the equations which describe a cir- 
cuit. Thus, computer programs for time domain 
circuit analysis have been among the important 
tools which have been developed to study and, 
ultimately, to control the effects of radiation 
on electronic circuits. CIRCUS was developed to 
meet such a need, and has been used extensively 
for the past three years in simulation, predic- 
tion, and radiation hardened design of complex 
circuits exposed to various radiation environ- 
ments . 


Although motivation for the development of 
CIRCUS originated from efforts to evaluate 
transient nuclear effects, the program is appli- 
cable to the much wider class of everyday 
electronic design problems, using a variety of 
different electrical forcing functions. In 
addition, since most physical systems have an 
electrical analogue, it may be extremely useful 
for investigating the transient behavior of 
mechanical, hydraulic, physiological, and other 
systems . 

This paper will present a brief description 
of some of the CIRCUS features of particular 
interest to the user. 

CIRCUS Features 

CIRCUS (Circuit Simulator) Is capable of 
analyzing any arbitrary circuit made up of 
resistors, capacitors, inductors, semiconductor 
devices, and voltage and current sources. As in 
other programs of this type, special efforts have 
been made in the design of CIRCUS to ensure ease 
of use and simplicity in the input required by 
the programmer. 

Like other time -domain circuit analysis 
programs such as NET-1, PREDICT, and SCEPTRE, 
CIRCUS requires a topological description of the 
circuit which specifies element values and 
forcing functions, and from this information 
constructs the appropriate time -domain equations, 
finds the de initial conditions, and then evalu- 
ates the transient solution by numerical integra- 
tion of the differential equations. Each program 
offers its special features, adaptability to 
given computer types and inherent advantages and 
disadvantages. It is assumed that the circuit 
designer wishing to use computer aids has in mind 
certain requirements for program features and 
flexibility when he chooses a given computer 
program for his specific needs. A complete 
description of CIRCUS is given in References 1 and 
2 and from a previous study was shown to be sig- 
nificantly faster than previous programs (Refer- 
ence 3)* In general CIRCUS performs the same 
type of circuit analysis as SCEPTRE, with special 
care taken to provide many user -oriented features. 
However, SCEPTRE has features, if needed, which 
make it a desirable tool for the development of 
equivalent circuits. A brief summary of some of 
the more important features of CIRCUS are the 
following: 
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X. Stored Semiconductor Models and Device Param - 
eter Library -- CIRCUS has stored device models 
for conventional and field-effect transistors ; 
conventional, tunnel, and Zener diodes; and 4- 
region devices. For each of these models, there 
are truilt-in phptocurrent generators for radiation 
effects work. Of course, these models, in combi- 
nation with other elements, can give rise to other 
equivalent circuits of interest. One has only to 
specify the terminal connections of a given device 
and either (a) request the program to obtain the 
device data from the Bevice parameter library 
tape, or (b) supply the appropriate device param- 
eters as part of the input, or some combination 
of both procedures. 

The stored models are based on a charge - 
control description of semiconductor device 
operation* These models in CIRCUS have been 
transformed into current rather than charge 
parameters and appear in an equivalent T-circuit 
very similar to the Ebers-Moll description. In 
fact, simple algebraic conversion can be made 
between charge -control parameters and Ebers-Moll 
parameters as described in (l), and this procedure 
has been used to demonstrate that identical cir- 
cuits analyzed on both NET-1 and CIRCUS with their 
respective stored models give an identical tran- 
sient analysis (Reference 3). Both these models 
ultimately reduce to the same differential equa- 
tions, and accurately describe both large and 
small signal operations in all four regions: 
cutoff, active -normal, active -inverted, and 
saturation. Storage effects are included. 

2. Special Purpose Parameter Variations and 
Multiple Analyses -- CIRCUS allows, in one com- 
puter run, multiple analyses of a given topology 
circuit with arbitrary changes in any or all 
element values, or analyses of many circuits with 
different topologies, CIRCUS also will allow a 
change In parameters (on a restart mode) of any 
given analysis, as well as at any predetermined 
time during an analysis. These latter capabi- 
lities have many uses, but among those most 
obvious are: (a) changes in parameters affecting 

short-time constants in low frequency circuitry 
after fast transients have passed, allowing for 
more rapid calculation; (b) degradation of device 
performance at a given time due to nuclear radia- 
tion; and (c) changes in forcing functions of 
circuits in AC equilibrium. 

3* Output Options — A variety of options are 
available as output for CIRCUS. In addition to 
any or all node voltages, and terminal currents 
for any elements, the user may request (l) in- 
ternal currents and junction voltages internal to 
any semiconductor device, (2) a monitor on voltage, 
current, or power dissipation in any element, 

(3) voltage across resistors, capacitors, 
inductors, and current sources, (4) the sum or 
difference of any two node voltages, or (5) any 
variable minus its DC steady state pulse. Only 
those variables selected by the user will be 
printed and/or plotted for display, and all in- 
formation is displayed at selected, variable 
intervals chosen by the user. These output 


options allow the designer to monitor any desired 
variable, and to perform optimization studies. 
They allow the circuit or system analyst to 
identify failure modes easily. 

4. Exponential Method of Numerical Integration - - 
CIRCUS solves the network differential equations 
for the transient analysis by using a variation 
of an exponential method developed by D. A. Pope 
(Reference 4). As programmed into CIRCUS, this 
integration routine is a one-pass, variable -step 
method controlled both by the solution behavior 
and by the smallest local time constant of the 
network. The fundamental step -size control 
parameter is related to the number of terms re- 
quired for convergence. Experience with a wide 
variety of circuits has shown that this exponen- 
tial method uses a step-size which is from three 
to four times as large as in the fourth -order 
variable -step Runge-Kutta method. The machine 
time per step is about the same in both methods 
since the time spent in the three extra deriv- 
ative evaluations per Runge-Kutta step is approx- 
imately the same as the time spent by the expo- 
nential method in computing the exponential of a 
matrix. Therefore, due to the larger step size, 
the exponential method will evaluate a given 
transient solution in one -fourth to one -third the 
time taken by a fourth -order Runge-Kutta method. 

5* Machine Independence -- CIRCUS has been coded 
almost entirely (98#)" 'in Fortran IV and thus can 
be used on any computer which has a Fortran IV 
compiler and adequate storage. CIRCUS has al- 
ready been coded for use on the following 
machines at several installations in the country: 
IBM 7090/7094, UNIVAC 1107/1100, G.E. 625/635, 
and the CDC 6600. It is expected that an IBM 
360 version will be available in March of 1967. 

Conclusion 

CIRCUS in its present form already repre- 
sents a flexible and efficient tool for computer- 
aided design. It is: 

1. user-oriented rather than programmer- 
oriented, meaning that its special 
features (multiple reruns with changed 
parameters, restart mode, stored models, 
etc. ) are designed so that the electronic 
engineer or vulnerability analyst will 
find it easy to use. 

2. computationally fast, and thus an 
economical tool as circuit analysis 
programs go. 

3* a proven and debugged program as 
attested by three years of un inter- 
rupted use. 

4. relatively machine -independent, and has 
already been coded for use on four dif- 
ferent computer systems. 

Ideas for further development of the program 
include new component models and more 
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computational flexibility. One specific compu- 
tational goal is to provide the capability for 
fast, efficient analysis of low frequency circuits 
having high frequency components. It is expected 
that resulting modifications to the program will 
significantly enhance the possibilities of 
computer-aided design, with the ultimate goal of 
increasing the productivity of the design engi- 
neer, and improving the performance and reliabil- 
ity of the final product at a significant reduc- 
tion in cost and time. 
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location on the circuit board, scaled' on - 
the 100 x 100 grid. Each element is a 
computer word and contains a number repre- 
senting the "altitude" of its location. 
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SUMMARY 

A new method has been programmed to 
produce the conductor path layout as part 
of the continuing development of the 
ACCEL system. The method is based on a 
topographic simulation scheme where com- 
ponent connector pin locations are repre- 
sented by peaks, and the areas between 
pins are represented by combinations of 
slopes, plains and valleys. Examination 
of this terrain can lead to the establish- 
ment of trails to represent conductor 
paths, Also, methods of displaying the 
simulated topography and trails have been 
programmed. Although this presentation 
is centered around the circuit board prob- 
lem, the technique may be useful in other 
problems of a contextual nature. 

Introduction 

The topographic simulation technique, 
presented here, is part of the continuing 
development of the ACCEL 1 system. It is 
based on the idea that the configurations 
of a circuit board can be represented in 
a topographic structure. This idea was 
first introduced to the authors by Dr. Iben 
Browning, Executive Director of the 
Thomas-Bede Foundation in Los Altos, Cal- 
ifornia, 

Topographic simulation, per se, is 
not new; It has been used in physical 
models to represent interrelated, steady- 
state forces. 2 However, it is believed 
that the sense in which it is used here is 
new; i.e., combining in the computer a 
topographic scheme with self -modifying 
characteristics in such a way as to pro- 
duce a dynamic topography. As used here, 
this dynamic property is especially vital 
in representing the constantly changing 
conditions caused by the growth of circuit 
paths. Another unique feature allows any 
given point to be "negative" in relation 
to all other points; i.e., everything is 
downhill from everything else. So, in 
essence, we have a multi-dimensional, dy- 
namic topography . 

Topographic Structure 

In the computer, the topographic area 
is represented by an array whose dimen- 
sions are 100 x 100 elements (this size 
is limited only by the available core 
storage in the computer) . The location 
of each element represents a corresponding 


The topographic structure that will 
be created within this framework repre- 
sents physical features occurring on the 
circuit board. At the beginning of the 
process, the following three features can 
be represented : 

1. The shape (outline) of the 
circuit board 

2. Holes or other obstacles to con- 
ductor paths 

3. Component and connector pin 
locations (lands) . 

The conductor paths, plated-through holes, 
and certain trouble spots that are estab- 
lished by the process are added to the 
topographic structure as they occur. 

The process begins with all elements 
being set at zero altitude. The data, 
representing electrical contact locations 
(component pins) and the locations of phy- 
sical obstacles (such as holes in the 
board) , are fed into the system. A "peak" 
is created in the grid at each such loca- 
tion. That is, the grid element repre- 
senting a pin location is given a high 
altitude, say 1000; the elements surround- 
ing it are given a lower altitude, say 750; 
and the next ring of elements still a 
lower altitude, etc. Then, to represent 
the edges of the circuit board, a "ridge" 
is created around the periphery of the 
grid. This ridge slopes inward from each 
of the four edges, forming a basin that 
contains the peaks. The actual form of 
the peaks and the basin is variable, being 
governed by control data used in the 
equation described below. 

Although there are several ways in 
which the structure may be numerically 
developed, the one used most extensively 
in the current study is of the following 
form: 

K 

A = — — — 

C x D + 1 
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where A is the altitude of any given 
element, 

C is a constant that controls the 
slope of the basin and peaks, 

D is the distance (on the grid) be- 
tween the element and the nearest 
edge of the grid or between the 
element:,, the center location of 
'the nearest peaks, and 

K is the constant that determines 
the magnitude of the altitude. 

Additional information can be repre- 
sented by replacing K with a function. 

For example, if a peak represents a high- 
voltage contact on the board, K could be 
a function of voltage and thus cause low- 
voltage paths to avoid that area. Or, 
similarly, if a location in the circuit is 
particularly sensitive to the capacitance 
effects of conductive paths, the K (for 
peaks in that vicinity) could be a function 
of sensitivity and thereby provide better 
dispersion of conductors. 

The variety of conditions that can be 
accommodated by this topographic structure 
makes it an especially pov/erful technique 
for the problem at hand and, perhaps, for 
other problems that involve the inter- 
relationships of many conditions. 

After the basin and peaks have been 
formed, the process of establishing cir- 
cuit paths begins. That process is dis- 
cussed in the following section of this 
report. However, changes in the topogra- 
phy occur during the process and should be 
mentioned here in order to complete the 
description of this subject. These changes 
do not occur in any particular sequence 
but, instead, occur as the requirement for 
them is determined by the system. The fol- 
lowing conditions cause the topography to 
be modified: 

1. As the locations of circuit paths are 
established, ridges are built along 
the paths. This, in effect, forms 
valleys between completed paths. The 
purpose is to help the uncompleted 
paths find their way around or be- 
tween completed paths. 

2. When a path runs into a trouble spot 
in the terrain, it tends to wander in 
a tight zigzag fashion, trying to 
find a way out. Trouble spots may 

be caused by conditions such as low 
areas surrounded by peaks, by very 
narrow valleys caused by ridges of 
several completed paths , or by a V- 
shaped or U-shaped bend in a comple- 
ted path. If a path gets trapped in 
such an' area, the altitude of that 
vicinity is progressively increased. 
This causes the path to modify its 


direction and, hopefully, find a way 

out. It also helps other uncomple- 
ted paths to avoid that area. 

Up to this point, the topography that 
has been described is somewhat comparable 
to natural topography. This comparison 
is at least close enough that features 
occurring in natural topography, such as 
ridges and valleys and peaks, are useful 
in describing the imaginary topography in 
this system. The one remaining topogra- 
phic characteristic, yet to be described, 
has no natural counterpart, and therefore, 
must be imagined in a different way. 

Perhaps the best way to view this 
feature is to imagine that the topography 
described thus far is built on a rubber 
sheet that is stretched on a frame. Then 
when a step is being taken (a step repre- 
senting an increment of a circuit path) 
from some location toward a target loca- 
tion, the rubber sheet is pressed down at 
this point of the target location, form- 
ing a cone with the ridges, peaks, etc., 
on its inside surface. When the step has 
been taken, its target location is re- 
leased and the system is prepared to find 
the next step. An iterative step-sequenc- 
ing system is used (that is, commutation 
of all of the uncompleted path ends) which 
gives the effect of developing all paths 
simultaneously. Therefore, the next step 
will probably be with an increment, a path, 
and a target different from those in the 
preceding step. The new target location 
is pressed down, a step is taken, and the 
target is released. This process of press- 
ing down and releasing successive target 
locations is an integral part of the pro- 
cedure for establishing circuit paths. 

Its exact purpose will become more evident 
in the following section, where the method 
of evaluating a path step is discussed. 

The following figures, produced by 
Fortran programs processed on the IBM 7094 
and plotted by the SC4020 plotter, help 
the reader to visualize the simulated 
topography. Figure 1 corresponds to a 
contour map showing the basin and peaks 
at the beginning of the process. Figures 
2 through 5 correspond to relief maps. 
Figure 2, like Figure 1, shows the basin 
and peaks; Figure 3 represents the same 
structure, except with one of the peaks 
pushed down to form a target cone. Fig- 
ure 4 contains the same set of peaks, but 
is later in the process after ridges have 
been built to represent completed paths. 
Figure 5 shows the effect of a target 
cone at this stage. 
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Algorithm for Path-Step Selection 

The system contains two general fea- 
tures that work together in the just- 
described topography to establish the lo- 
cation of Circuit paths. One is the basic 
algorithm for step sequencing and step 
sequencing and step evaluation described 
below. The other comprises the usage 
and variations of this algorithm in the 
several phases of the system which are 
described in the following section. 

The definition of a path step, as 
used here, is an element in the 100 x 100 
grid and represents a portion of a con- 
ductor path. If a path is to be developed 
between peaks A and B, it will be composed 
of at least as many steps as there are 
elements in the grid between A and B. 

Since paths frequently need to wander 
(for example, around other peaks), the 
number of steps to complete a path de- 
pends upon the obstacles between A and B. 

The step-selection procedure is 
basically one of examining the slope of 
the terrain in the immediate vicinity of 
the element from which a step is to be 
taken and choosing the direction that is 
most downhill. The exceptions to this 
rule are that a step must not be on or 
adjacent to steps of a different node, 
and a step cannot go backward (i.e., a 
path cannot intersect itself) . Because 
of these exceptions, the most downhill 
direction may be unacceptable. Therefore, 
it is necessary to provide alternate 
choices of step direction. Another 


restriction, which is purely for the sake 
of simplicity, is that a step must be one 
of the four orthogonal elements of the 
current location (i.e., it cannot step 
to a diagonal element) . 

An interesting feature of the step- 
ping algorithm is that the two active ends 
of a path try to find each other. This 
stepping is controlled by a data array 
that keeps track of the most recent step 
location of each path end, as well as the 
point of origin (peak) of each path end. 

For definition, the most recent step of 
a path end is called the "target" of the 
other end: the peak toward which a path 

is moving is called the "target peak" 

(i.e., the point of origin of the target 
path end) ; and the path end that is in the 
process of taking a step is called the 
"active end." During the step-evaluation 
process, three characteristics of altitude 
are measured for each of the orthogonal 
elements at the active end: 

1. The actual altitude currently in 
the topography. 

2. The portion of item 1 which is 
due to the target peak, visual- 
ized as though the target peak 
had a cumulative effect upon 
the entire terrain. 

3. The effect on altitude caused by 
the target being pulled down to 
form a cone, as described earlier. 
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Items 2 and 3 are computed values pro- 
duced by the following equations: 

Altitude due to Target Peak = — ■ — 

C x D + 1 

where C is a constant that defines the 
shape of the peak, and D is the distance 
(grid elements) between the center of the 
peak and the element being considered. 

Altitude Change due to Cone = 

K x (100 - D) where K is a constant 
that defines the depth of the cone, 
and 

D is the distance from the 
target to the element being con- 
sidered. 

Then, for each of the orthogonal elements, 
an effective altitude is computed as be- 
ing item 1 minus item 2 minus item 3. 
Whichever of these elements has the low- 
est effective altitude and does not vio- 
late any of the restrictions described 
earlier will be the new step location. 

Process Phases -- The overall operation 
of the system is divided into several 
phases. Each phase, while employing the 
algorithm previously discussed, follows 
a slightly different policy in making 
connections. At the present time, the 
system employs four separate phases. A 
phase control subroutine is used to over- 
see the sequence of phases. It is easily 
varied to try different policies for mak- 
ing connections. 

Front of Board — Phase 1 — Phase 1 begins 
with a topography that includes only the 
border and the peaks at the land locations. 
During this phase, steps are taken from 
both ends of a connection pair toward 
each other. A step from both ends of all 
connection pairs constitutes one cycle of 
the stepping algorithm. It is important 
to note that during this phase there are 
no ridges to represent paths in the ter- 
rain. However, a step is prohibited if 
its location would be adjacent to a pre- 
vious step from a wire on a separate node. 
Thus there is an infinite barrier sur- 
rounding the wires as they progress, even 
though the sloping ridges are not present. 
Consequently, no early warning is provided 
to warn of an impending obstruction. Af- 
ter the cycle limit is reached or all the 
wires become trapped, whichever comes 
first. Phase 1 terminates. Typically, 
about half the connections are made dur- 
ing Phase 1. 

Phase 2 — Phase 2 is almost a repeat of 
Phase 1, but with one important differ- 
ence. All the wires completed in Phase 1 
are provided With ridges, sloping down- 
ward and outward from both sides of the 


connection path. Of course, all those 
connections not completed are wiped out, 
and Phase 2 attempts them once more, but 
now with more contextual information from 
which to work. 

Phase 2 is the last attempt to make 
the connections on the front of the board. 
Therefore, during this phase, a special 
routine is used to determine whether the 
steps from the two ends of the wires are 
still making progress toward each other. 

If this routine determines that one end 
is beginning to stray, its progress is 
halted and a land is placed there. Thus, 
when Phase 2 terminates, for each wire 
that has still not completed a connection, 
a land is placed at the last position of 
the two ends. 

If a two-sided board is being gener- 
ated, the program now tries to complete 
the remaining connections, between the 
new land locations, on the back of the 
board. The completed path locations that 
have been established for the front side 
of the board are saved and the topography 
is initialized to zero before beginning 
on the back side. 

Back of the Board — Phases 1 and 2 — The 
topography representing the back of the 
board is developed, placing the border as 
before, and creating the peaks as before. 
However, there will be many more lands or 
peaks than there were for the front of 
the board. The newly generated lands at 
the ends of uncompleted wires, as well as 
the original lands, make up the data for 
calculating peaks. (These land positions 
must be avoided because the lands pro- 
trude through the circuit board.) 

After setting up the terrain for the 
back of the board. Phase 1 is utilized as 
before. Of course, only those wires still 
uncompleted are attempted , and they step 
from new land locations. 

Again, after Phase 1 is completed, 
dikes are generated for newly completed 
paths, and uncompleted paths are erased. 
Phase 2 is then tried. After Phase 2, if 
there are still uncompleted paths. Phases 
3 and 4 are tried. 

Phases 3 and 4 — Phases 3 and 4 differ 
from Phases 1 and 2 in two ways: (1) only 
one wire is tried at a time and (2) steps 
are taken from only one end of the wire. 
Thus, these phases are executed once for 
each uncompleted wire. Phase, J3 is tried 
first. During this phase, end 1 of the 
wire is trapped. In other words, it stays 
at its land location while only end 2 
steps toward that opposite land. It is 
important to note that during this phase 
steps are taken for only one wire path. 
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Figure 6. Stepping Sequence, Phase 1 


Figure 7. Stepping Sequence, Phase 2 


After all wires have been tried. Phase 4 
is entered. It is identical to Phase 3, 
except that the remaining uncompleted 
wires are tried from the other end. That 
is, end 2 is trapped, while end 1 does the 
stepping. After all wires that are still 
uncompleted have been attempted by Phase 
4, the program generates a plot tape for 
producing a picture of the board and con- 
ductor paths and then goes on to the next 
circuit to be routed. 

The following figures illustrate 
the progression of paths in the terrain 
shown in Figures 1 through 5. Figures 6 
and 7 show the progress of paths in Phases 
1 and 2, respectively. Note that the first 
frame of Phase 2 contains those paths that 
were completed in Phase 1. 

The following figures illustrate the 
final product of the system, using both 
sides of the board. Figure 8 is a diagram 
of a particular circuit, showing the place- 
ment of components and the connections 
that need to be made. Figures 9 and 10 
are plots, produced by the system, to pro- 
vide the required connector paths. Fig- 
ure 9 represents the front side of the 
board, and Figure 10 the back. In this 
particular case, all connections are com- 
pleted. 

The component placement shown in 
Figure 8 was computer generated, using the 
current ACCEL system. 



Figure 8. Diagram of Circuit, Showing Placement of Components 
and of Connections to be Made 
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Conclusions — Although it is somewhat 
premature to compare the performance of 
this technique to the routing technique 
used in the operating ACCEL system, com- 
parative results from two test circuits 
reveal that, in terms of the number of 
connections completed, the topographic 
system did slightly better in both tests. 
Also, the arrangement and dispersion of 
paths appear to be better. On the other 
hand, the topographic system left a few 
paths uncompleted that the ACCEL system 
would have completed with ease. There- 
fore, a combination of the two techniques 
seems most promising. Since electrical 
characteristics (aside from connection 
requirements) can be represented in the 
topography, the system has distinct ap- 
peal over analytical methods that only 
consider minimizing crossovers. 
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SUMMARY 

This paper presents a balanced discussion of 
SCEPTRE (System for Circuit Evaluation and Predic- 
tion of Transient Radiation Effects) from both the aca- 
demic and practical points of view* Two basic inci- 
dence matrices that apply to any connected electrical 
network composed of resistances, capacitances, in- 
ductances, voltage sources and current sources will be 
used to derive a general F matrix. The F matrix will 
be used in turn to establish very practical relations 
between the tree branches and links of the general elec- 
trical network. It will be shown that this relationship 
is exploited to form an efficient algorithm that iB used 
in forming some of the equations that are required for 
a general transient solution. 

The second part of this paper will describe the 
most important features that are available to aid in the 
practical solution of contemporary networks. Each of 
these features were designed to provide more conve- 
nience to the user without compromising the overall 
flexibility of the program. 

INTRODUCTION 

Automatic circuit analysis programs are intended 
to relieve the user of the tedious and error-prone tasks 
of writing and programming circuit equations. As a 
result, the user is only required to enter a topological 
and quantitative description of each network of interest, 
according to an easily learned format. The general ac- 
ceptance of programs of this type can be gauged by the 
steadily increasing volume of reports and descriptions 
that are currently appearing in many technical journals. 

In 1963, a group at IBM's Electronics Systems 
Center was funded by the Air Force Weapons Labora- 
tory at Kirtland Air Force Base, New Mexico, to 
^The work described in this paper was spon- 
sored by^the Air Force Weapons Laboratory 
unde s^Sbntr act AF29{60 l)-b852/ 'Supplemental 
Agreement 1, as part of the Defense Atomic 
Support Agency program under NWER, Subtask 
16. 015. 


develop an automatic circuit analysis program. This 
program waB to be flexible enough to accommodate a 
variety of hoh-conventional effects that are of interest 
in transient radiation problems as well as general 
enough to handle conventional transient problems. The 
result was PREDICT (Prediction of Radiation Effects 
Using Digital Computer Techniques), an automatic 
transient program that has been in general dissemina- 
tion throughout the country since 1964. Inevitably, 
the heed for a successor program (SCEPTRE) slowly 
became evident. Invaluable comments, pro and con, 
flowed in from the more active and experienced users. 
This information played an important part in the de- 
velopment of many of the features that today distin- 
guish SCEPTRE from PREDICT. The reader who is 
seriously interested in the progress that has been made 
in this field would do well to read and compare the re- 
spective manuals for these programs (PREDICT [l] , 
SCEPTRE [2] ), Dissemination of both programs is 
controlled by Lt. Gary Pritchard, Air Force Weapons 
Laboratory, Kirtland AFB, New Mexico 87117. 

Symbols 

Y vector of state variables 

Y vector of state variable derivatives 

b number of network elements 5 i. e. , sum of 
all resistors, capacitors, inductors, 
voltage sources and current sources 

n number of network nodes 

1^ vector of currents through all network 
elements 

Vjj vector of voltages across all network 
elements 




V 
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F an incidence matrix that expresses a rela- 
tion between network links and tree 
branches 

T 

F the transpose of F 

The description of a SCEPTRE formulation, must 
first begin with the definition of a tree [3]. A tree 
is a connected subgraph of a network which includes 
every node in that network, but contains iio closed loops 
dr meshes. A specific tree may be defined as one in 
Which elements are inserted in the following order of 
priority: voltage sources, capacitors, resistors, and 
inductors. All current sources are excluded from this 
tree. All elements that are included in the tree are 
called tree branches; all elements that are excluded 
from the tree are called links. 

The state variables of any system may be defined 
as the minimum set of variables which, together With 
all inputs , are sufficient to determine all other system 
quantities at any particular instant of time. All system 
quantities at time tj, are dependent upon the values of 
the state variables and inputs that hold at time tn and 
are not directly dependant on those at any other time. 

It can be shown that the state variables of the general 
electrical network are the set of capacitor tree branch 
voltages and inductor link currents. This choice of 
state variables is equivalent to specifying the independ- 
ent initial conditions of the network. It may be stated 
further, that the choice of state variables, leads to a 
set of simultaneous first-order differential equations 

M • 

The solution process is best illustrated with the 
highly simplified block diagram of Figure 1. The pro- 
cess begins with the insertion of the state variables Y 
that are valid at time where n = 0 at the start of a 
transient problem. 1 All currents, Voltages and state 
variable derivatives that are valid at t =t h are computed 
and output (if desired). The state variable derivatives 
Y (tn) are numerically integrated by the integration 
routine to provide the state variables Y (t n +j) that are 
valid at the next solution increment. It is worth noting 
that the block diagram serves to illustrate the two most 
basic factors concerned with the amount of computer 
time required for the solution of transient problems. 
The number of solution steps required to complete a 
given problem depends on the integration routine used, 
and the amount of computer time required per solution 
step depends on the efficiency with which the algebraic 


‘ The initial values of the state variables may be 
supplied, if known, by the user. Otherwise, they may 
be automatically determined by an iterative procedure 
that is contained in SCEPTRE. 


operations are carried out in the Equation Solution 
block. 



Figure i. Solution Process Diagram 

The topological makeup of any connected network 
may be described by means of two incidence matrices. 
The first of these is the fundamental cut set matrix Q 
where: 

Q = [qij] is a matrix containing n-1 rows and b 
columns for the general connected network con- 
sisting of n nodes and b elements where: 

q t j = +1 if the direction of fundamental cut set 
i coincides with the reference direction 
of element j in fundamental cut- set 1 

q^j = — 1 if the direction of fundamental cut-set 
i is in opposition with the reference dir- 
ection of element j in fundamental out- 
set i 

q^j = 0 if the fundamental cut-set i does not 
include element j 

The direction of each fundamental cut set Is custom- 
arily chosen to agree with that of its associated tree 
branch. It will be convenient to partition Q according 
to links and tree branches as 

(1) [QL Qt] 

which, if the rows and columns are ordered appropri- 
ately, becomes 

Q - [“F t tl] 

where the (n-1) rows and columns of the unit matrix U 
correspond to the tree branches and the (b-n+1) columns 
of -F t correspond to the links associated with any tree. 

Another incidence matrix can be defined which rep- 
resents the topological relationship between the ele- 
ments and the meshes, or loops, or circuits, of the 
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network. This matrix will be called the fundamental 
circuit matrix T, where 

T = ftijl is a matrix containing (b-n+1) rows and 
B columns for the general connected network 
consisting of n nodes and b elemehts where: 


is always true, a simple partition permits 
(8) [u F] 

or 



tjj * +1 if element j is part of fundamental 
circuit i, and if the reference directions 
of element j and fundamental circuit i 
coincide. 

ty = -1 if element j is part of fundamental cir- 
cuit i, and if the reference directions of 
element j and fundamental circuit i do 
not coincide 

tij = O if element j is not part of fundamental 
circuit i 

The direction of each fundamental circuit is custom- 
arily chosen to agree with that of its associated link. 
The number of rows in the T matrix corresponds to the 
fact that only (b-n+1) circuit equations are required for 
the solution of any network. If the columns of T are 
ordered appropriately, T may be partitioned as 

(2) T = [u f] 


(9) V L = -F V T 

where, V, and V T are the vectors of link voltages and 
tree branch voltages, respectively. Equation 9 is 
significant because it isolates the link voltages of the 
general network in terms of the tree branch voltages. 
The Q and T incidence matrices have been eliminated. 

Equations 6 and 9 indicate that the F matrix can be 
used to transform the vector of link currents into tree 
branch currents and the vector of tree branch voltages 
into link voltages. These relationships also indicate 
that F contains one row for each network link and one 
column for each network tree branch. If the tree of the 
sample network of Figure 2 is indicated by the bold 
lines, the appropriate F matrix would appear as in 
Figure 3. The associated elements are used to caption 
each row and column. Note that the three link voltages 
may be "read” horizontally from the F matrix as: 

VR1 = -VC1 + El 


where, the (b-n+1) columns of the unit matrix U cor- 
respond to the links and the n-1 columns of F corres- 
pond to the tree branches of any tree. If the columns 
of Q and T are arranged in the same element order, it 
is well known that: 

(3) Q T T = 0 


Also, since 

(4) Q I b = 0 


it always true, a simple partition permits 


(5) [-F T U] 



= 0 


VL1 = VC1 - VR2 
VJ1 = VR2 


R1 LI 



leading to 

(6) I T = F t I I 

where II and It are the vectors of link currents and 
tree branch currents, respectively. Equation 6 is 
significant because it isolates the tree branch currents 
of the general network in terms of the link currents. 

In the same way, since 

(7) T V b = 0 


TREE BRANCHES 


LINKS 



Figure 3. F Matrix for Sample Network 
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These equations are just a restatement of matrix equa- 
tion 9. Note also that the three tree branbh currents 
may be "read" vertically from the F matrix as: 

IC1 = IR1 - IL1 

IR2 = IL1-J1 

IE1 = IR1 

These equations are a restatement of matrix equation 6. 
SCEPTRE is programmed to "read” the F matrix that 
is constructed for each network under analysis and to 
store the resulting equations^ for the duration of the 
problem. These stored equations may then be updated 
at each solution step. This method may be contrasted 
with the more conventional approach used in older pro- 
grams (PREDICT, for example) in which, the same 
results are obtained by repeated manipulation of the 
matrix equations at each solution step. When one con- 
siders that many problems require literally thousands 
of solution steps, the computational savings that may 
be gained by this process are quite significant. 

Some of the features of SCEPTRE that represent a 
significant advance over PREDICT will be found to be 
unique, or nearly so, if one canvasses all the automatic 
programs that are available today. It should be real- 
ized that this discussion serves only as a summary; a 
much more complete description exists in [ 2 ] . 

Optional Stored Models 

Most contemporary circuits include various types 
of diodes and transistors that require equivalent cir- 
cuits with varying degrees of complexity. For example, 
a large signal equivalent circuit for a transistor usu- 
ally requires four to six passive elements and at least 
three to four current generators. If the computer pro- 
gram lacks a stored model feature, the user has no re- 
course but to insert each component individually into 
the overall circuit. A partial remedy for this problem 
is the provision of an "exclusive" stored model of con- 
ventional form for diodes and transistors. The word 
"exclusive" is used in the sense that the user must use 
whatever topology is originally provided in the stored 
model; no provision is made to allow topological choice. 
A far better solution is called the optimal, stored model 
feature. The user may himself, by a simple procedure, 
store almost any3 equivalent circuit that is composed 


2 

’The topological limitations that apply to this method 
are given in [ 5 ] . 

3 

* Network configurations that include voltage source 
loops or current source cut sets are prohibited. 


of; resistors, capacitors, inductors (including mutual) 
inductances), and voltage and current sources. Once 
stored, any model may be called out for use as many 
times are required in a given circuit application. This 
feature is not limited to equivalent circuits of active 
devices; frequently used subnetworks such as biasing 
networks and filter sections may also be stored. 

Rerun Capability 

It is often desirable to know how well a given net- 
work would function if one or more of its components 
or forcing functions were varied in size. Separate runs 
can always be made to determine this, but a better way 
now exists. A master run is first described, which 
includes the nominal values of all components and forc- 
ing functions. The user appends to this description a 
list containing the number of reruns desired and the 
elements that are to be changed for each rerun. No 
topological changes are permitted between the master 
run and any of its associated reruns. 

Since certain operations are common 4 between the 
master run and its reruns, machine time economy can 
be effected by simply performing these operations once, 
for all of the runs. It should be mentioned, however, 
that this rerun facility does not imply that anything like 
true statistical analysis is economically feasible. So- 
lution time per run is still often measured in minutes 
and therefore, a large number of reruns can become 
very expensive. 

Automatic Termination 

In some instances, the user will have no interest 
in carrying a transient run to completion if it is known 
that some voltage or current has exceeded, equaled, 
or dropped below some other network quantity or cri- 
terion. A single entry into the program can set up the 
logic necessary to monitor and automatically terminate 
any run. Practical applications could be concerned 
with turn-on, turn-off, or saturation of a transistor or 
any specific current or voltage limit that the user would 
care to set. This feature may be used concurrent with 
the Optional Stored Model and Rerun features. 

Flexibility 

Some idea of the trend toward flexibility has already 
been discussed in the section on Optional Stored Models, 
but there is considerably more to the subject. Sources 
may be inserted at almost any point in a network and 
multiple independent variables may be used to control 


4. 

The tree structure, F matrix, and basic equation 
setup do not change between runs unless the network 
topology is changed. 
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these sources. In addition, a special section has been 
built into the program to permit the definition of pseudo 
circuit parameters. To illustrate, the user may deter- 
mine the power dissipation in resistor R1 of any network 
by defining the parameter PR1 = (R1 * IR1 ** 2). The 
indicated mathematical operation would be performed 
and output at each solution increment. The same fea- 
ture allows the user to solve systems of first order dif- 
ferential equations that may or may not have any connec- 
tion with an electrical network. Consider the system 

W1 = 3W1 + 2W2 + 7 

W2 « W1 - 3W2 + 3W3 

W3 = 10W2-20W3 + 5 

where, Wl(t Q ), W2(t 0 ), W3(t Q ) are all known. The 
complete transient response to this system can easily 
be obtained with just a few input data cards. 

Subprogram Capability 

It appears to be obvious that no program can be de- 
vised to accommodate every operation that all users 
may require. The next best solution is to write the 
program so that users can insert subprograms of their 
own design which will perform special purpose tasks. 
This option is available to SCEPTRE users with the 
stipulation that all such subprograms are subject to the 
rules and requirements of Fortran IV. 

The preceding paragraphs by no means exhaust the 
number of special features that have been incorporated 
into SCEPTRE, however they should serve to convey 
some idea of the direction in which automatic circuit 
analysis is moving. From all indications, development 
of programs of this type will continue at a rapid rate 
in the foreseeable future. 
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SUMMARY 


Accurate and timely circuit engineering is based upon: (a) accurate circuit parameter characterization* 

(b) flexible network analysis computer programs, (c) continuous computer program development, and (d) empirical 
investigations in adequate laboratories. This paper shows how these criteria are applied to obtain meaningful cir- 
cuit analysis results. Typical examples of circuit performance obtained through computer -assisted analyses are 
given and compared to measured results . 


INTRODUCTION 

This paper describes the capabilities an engineer- 
ing organization must consider to perform circuit de- 
sign and analysis efficiently with the aid of computers. 
The availability of a large computer system is 
assumed . 

Circuit engineering considers the worst-case cir- 
cuit performance , which is a function of piece-part 
tolerances, temperature environment, and aging. Cir- 
cuit performance criteria include dc stability, power 
dissipation, and dynamic performance. Dynamic per- 
formance encompasses transient behavior, dynamic 
stability, and four -terminal network characterization 
of module circuitry. While emphasis will be given to 
circuit analysis and design, system design and analysis 
can be and has been performed using the basic 
approaches described herein. 

Accurate and timely circuit engineering is based 
upon: (a) accurate circuit parameter characterization, 
(b) flexible network analysis computer programs , (c) 
continuous computer program development, and (d) 
empirical investigations in adequate laboratories. 

ACCURATE CIRCUIT 
PARAMETER CHARACTERIZATION 

Overall accuracy in the computer programs used 
emphasizes the need for accurate parameter data of all 
piece parts , information that is seldom available from 
device manufacturer's specifications. Models giving 
accurate device representation can be used because of 
the increased capability afforded to the engineer by the 
computer. Device parameters , therefore, must be 
chosen to fit the models and then measured. For tran- 
sistors and diodes these parameters are measured on a 
limited sample (25 devices) at 25 degrees centigrade, 


and statistical techniques are used to approximate 
worst -case parameter limits . Data extrapolation for 
other operating temperatures is accomplished by semi- 
empirical relationships , developed to reduce parameter 
measuring costs . Integrated circuit parameters are 
measured, depending on the module configuration, and 
are defined as four -terminal network parameters on 
analog devices, or as switching times and amplitudes 
for digital devices. In special cases, integrated circuit 
parameters can be measured by gaining access into the 
chip; but, after the measurements, the device is no 
longer usable . 

For dc and transient analysis , the nonlinear tran- 
sistor model (Figure 1) that is based on the Ebers-Moll 
large signal switching model is used. The model is 
valid for the three regions of transistor operation: 
cutoff, active, and saturation. The parameters re- \ 
quired, I E vs V BE , I CB vs V BC , C b , 0 and C b , e , are 
readily obtained through laboratory measurements. 

The hybrid-Pi ac small signal transistor model 
shown in Figure 2 indicates the circuit configuration, 
model parameters, and the semi -empirical equations 
relating temperature dependence . Since the model para- 
meters are interdependent, they must be defined in 
terms of the basic independent parameters to determine 
worst -case parameter variations . The basic parameter 
variations are obtained from the common -base and 
common emitter h parameters measurements , which 
are accurate and readily obtained within the laboratory. 

NETWORK ANALYSIS COMPUTER PROGRAMS 

Experience has shown that circuit analysis pro- 
grams, to be useful in circuit engineering, must have 
the capacity to handle large and complex circuits , such 
as the circuits shown in Figures 3 and 10. Four basic 
programs are used, depending on the circuit 
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Table I —Program Capabilities 



Base - Emitter Voltage 




AC-DC 

PANE/360 

TRANSIENT 

CONTROL 

SYSTEM 

S-PLANE/360 

ANALOG 

SIMULATOR 

DSL-90 

INPUT 

PARAMETERS 

200 

200 



OUTPUTS 

50 

100 


50 

DEPENDENT 

NODES 

60 

60 

20 


DISTRIBUTIONS 

SHAPES 

6 

__ 




TABLES FOR 
NONLINEAR 
INPUTS 

20 

15 



TRANSFER 

FUNCTION 

ORDER 



25 


FORCING 

FUNCTION 

— 

__ 

5 

__ 

STANDARD 

FUNCTION 

GENERATORS 




25 

SIMULATION 

BLOCKS 

— 

... 


50 

MEMORY 

REQUIREMENTS 

(BYTES) 

256K 

220K 

200K 

150K 


•PREDI CT/300 incorporates all l So problom-iolvln* capabilities of PREDICT I, which w«» developed by IBM 
uodor contract AF!9(tioi)-6399 for the Roaoaroh and Technology Division. Air Force Systems Command, 

Air Fore* Woapons Laboratory, Kirkland Air Fores Base, Now Mexico. 


c 


E 


C 



Figure 1« -Nonlinear Transistor Characteristics and Equi- 
valent Circuit, DC and Transient Analyses 

performance required. The program characteristics 
are given in Table I. 

For linear circuits , the Performance Analysis of 
Electrical Circuits program (PANE) is used to obtain 
dc stability , bias levels , peak and average power dissi- 
pation, and frequency and phase response. The circuit 


analysis is performed automatically by the program 
from a topological circuit description, requring a mini- 
mum of programming knowledge. Tabular and graphi- 
cal outputs of worst -case performance or statistical 
performance are obtained. "Piece-wise" linear dc 
analysis is performed, using this program, for analysis 
of switching circuits to determine ON and OFF states , 
and peak and average power dissipations. 

PREDICT/360, a transient analysis program, is 
capable of computing circuit responses , as a function 
of time, to any periodic or nonperiodic forcing function 
which can be defined by a table or equation. The pro- 
gram outputs are time history tabulations or plots or 
both. The program in its present form is not a worst- 
case analysis program and engineering judgment has 
been used to determine circuit worst-case performance . 
By considering the worst- case parameters that would 
degrade or improve the circuit and by rerunning the pro- 
gram, the worst-case performance, limited only by the 
degree of engineering judgment, can be determined. 

Linear circuit or control system dynamic response, 
including root -locus analysis , is obtained by the S -PLANE 
program. Similar to PREDICT, S -PLANE obtains the 
transient response of a transfer functicm to an impulse , 
sine wave, square wave, ramp, or exponential input. 

The transfer function of a network can also be obtained. 
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Semi-empirical Equations Showing 
Temperature Dependence 
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Note: "T" is in degrees centigrade . 


Figure 2. -*Hybr id- Pi Equivalent Circuit 

*Bymers R. J. "Linear Amplifier Design Manual," M05- 
001 1 -0, IBM, Standards Engineering, General Products 
Division, Glendale Laboratories 


The dc amplifier (Figure 3) is an analysis example 
demonstrating the performance characteristics deter- 
mined, using computer-aided circuit analysis techni- 
ques . The analysis determined the amplifier open- 
loop performance characteristics and, from these 
characteristics, a mathematical model was developed. 
The model was then used to determine the amplifier 
closed-loop performance . The worst-case dc null off- 
set was determined for the closed-loop amplifier as a 
function of temperature . Figures 4 and 5 show the 
computed amplifier open-loop voltage gain magnitude 
and phase along with laboratory measurements for 
comparison of computed and measured results . The 
amplifier worst -case open-loop voltage gain magnitude 
response is shown in Figure 6. Results of the worst- 
case analysis of the amplifier indicated a marginal de- 
sign with respect to amplifier stability. The amplifier 
was then stabilized with additional compensation. The 
characteristics of the improved amplifier are shown in 
Figure 6, depicting a worst -case phase margin of 40 
degrees. Figure 6 also shows the results of the dc null 
offset analysis and the open-loop-mathematical model 
derived. The model performance characteristic is also 
plotted for comparison. The closed-loop performance 
of the amplifier is shown in Figures 7,8, and 9. The 
circuit configuration analyzed is shown in each figure , 
as is the minimum, nominal, and maximum circuit 
response . 

For circuits or systems having highly nonlinear 
elements (such as level detectors in combination with 
amplifiers) , the circuit is broken into individual blocks 
that can be analyzed individually by the programs 
mentioned. It is this division that requires a high de- 
gree of engineering skill, judgment, and experience on 
the part of the analyst. These blocks are then charac- 
terized as four -terminal networks including all non- 
linearities, thus allowing a system simulation to be 
obtained by using the analog computer simulation pro- 
gram DSL-90. 

COMPUTER PROGRAM DEVELOPMENT 

While some of the programs stress that the lack of 
programming knowledge does not impede their use, 
experience has shown that being able to modify the pro- 
grams to fit a particular circuit under consideration has 
produced timely and quite unexpected results . Because 
of the increased flexibility afforded by the programs , 
the engineer can evaluate his design or complete his 
analysis, without ignoring important but often neglected 
parameters, if he can fit the program to the circuit 
under scrutiny. This program knowledge and the use of 
these programs have resulted in a continued program 
updating to enable the analysis programs to keep abreast 
of circuit technology. In addition, this knowledge has 
resulted in continued work to reduce computing time 
and to increase the flexibility of the programs (i.e. , 
adapting steepest descent techniques to PREDICT/360 
now under consideration, to enable worst-case transient 
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Figure 4. -Gain of Open-Loop Amplifier 
Magnitude (-55°C) 

analysis)., Program development over the past four 
years has increased from the analysis capability of cir- 
cuits having a maximum of 11 dependent nodes on a sta- 
tistical basis to the capability- afforded by the PANE pro- 
gram. Along with the ac and dc capability of PANE, 
other programs were introduced and updated to their 
present capabilities as shown in Table I. 



c .01 .1 1 10 100 IK 10K 100K 1M 10M 100M 

* FREQUENCY IN HERTZ 


Figure 5. -Gain of Open- Loop Amplifier 
Phase Angle (-55°C) 

EMPIRICAL INVESTIGATIONS 

Often, especially during design efforts when timely 
analysis results are required, empirical circuit re- 
sponses must be obtained, since mathematical models 
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Figure 3 1 . -Descriptive Block Diagram of 50 MA Servo Amplifier 


characterized using empirical techniques . The empir- 
ically determined load characteristic, driving point 
impedance function, and synthesized circuit used for 


the analysis are shown in Figure 12, The magnetic 
amplifier input impedance is shown in Figure 13. After 
| all blocks were characterized, the amplifier was 
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analyzed for closed-loop performance. From the 
closed -loop performance characteristics , a mathema- 
tical model for die amplifier was derived as shown in 
Figure 14 . The worst-case current gain performance 
characteristics are shown. The broken line depicts 
the minimum gain at the temperatime design limit of 
the amplifier. Laboratory measurements are plotted 
to show correlation between computed and measured 
results . 


CONCLUSIONS 

Computer-assisted circuit engineering, therefore, 
requires more than just one computer circuit analysis 
program. The circuit parameters must be known. 
Several programs must be available and checked out in 
the computer laboratory. Program development must 
keep the programs compatible with circuit technology . 
Finally, a laboratory for parameter measurement, 
empirical evaluation of circuit blocks , and anatysis 
reasonableness checks must be an integral part of 
computer circuit engineering. 
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Figure 12.-50 MA Servo Amplifier Load 



Figure 13. -Magnitude of Input Impedance of 450 
Turn Control Windings 
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Figure 14.-50 MA Servo Amplifier WC Closed-Loop 
Performance 
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SUMMARY 

The widespread availability of high- 
speed digital computers has profoundly 
affected many technical areas, and the 
impending introduction of powerful time- 
sharing systems will again revolutionize 
the engineer's work. One field of con- 
siderable importance is the area of 
optimization or approximation; for 
example, design of loss, phase, and delay 
equalizers is usually accomplished by 
iterative techniques. Many optimization 
algorithms have been designed and tested. 
Each possesses particular advantages and 
disadvantages. Search techniques such 
as grid, random, and pattern search and 
the simplex algorithm are less affected 
by certain function anomalies than slope- 
following methods such as steepest 
descent, gradient par tan, Fletcher -Powell 
and generalized Newton-Raphson. Steepest 
descent and gradient partan usually work 
with a poor initial estimate of the 
solution, while the generalized Newton- 
Raphson approach is efficient near the 
optimum, but may diverge elsewhere. The 
Fletcher -Powell technique utilizes both 
first and second derivative information, 
and is quite often the most efficient. 

Not all of these methods allow ready 
inclusion of constraints. Ideally the 
designer should have available a variety 
of algorithms, and should be able to 
apply any combination to his problems. 

INTRODUCTION 

The digital computer makes practical 
iterative design and optimization proce- 
dures where only an engineer's judgment 
could have been used previously. Although 
overall system optimization is not yet 
achievable, subsystem optimization is 
used every day in practical engineering 
problems. The introduction of direct- 
access time-sharing computer systems will 
again increase the range of problems 


wherein iterative optimization techniques 
are useful. 

In approaching an iterative approx- 
imation or optimization problem the 
engineer has many techniques from which 
to choose, varying from some brute force 
search techniques to highly sophisticated 
second order gradient algorithms. Each 
approach has its advantages and disadvan- 
tages. Several techniques will be briefly 
described in this paper, followed by a 
description of an existing program with 
sample results. 

For more detailed descriptions and 
comparisons of optimization algorithms 
the reader is referred to a number of 
excellent articles^- - ^ and books7-9 t 
Reference 8 presents results of recent 
practical work, comparisons of methods, 
and an unusually complete bibliography. 

The contents of this paper have been 
presented in expanded form. 10 

THE OPTIMIZATION PROBLEM 

The function to be optimized is of 
the form 

y - y(f 1> ...,f_ t , x n j (l) 

where f_ is the set of functional param- 
eters and x is the set of parameters 
which may be varied to optimize the 
function y. Usually y is to be maximized 
or minimized; in this paper it is assumed 
that y is to be minimized. 

In general y can be any function of 
£_ and x that is to be minimized. In 
practical optimization problems y usually 
takes the form of an error measure which 
is a complicated function of £_ and x. 

One popular functional form is the so- 
called "least squares" error criterion 

y = y w i (g i -r i ) 2 = y v ± e ± 2 (2) 

t r 
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where 

y = error measure to be minimized, 
e^ = unweighted error at £., 

*= g(f. ,x^, . . . ,x^) = computed values, 

x. = j th parameter of the n parameters 
^ to be varied, 

r. - r(f .) - required or desired 
function value, 

f . = i th value of independent 
variable, and 

w. = assumed weighting function 
(allows placing emphasis on 
particular requirements) . 

Another useful criterion is the maximum 
absolute error (Chebyshev error) : 

y = max | . (3) 

In both cases the function y is often a 
complicated expression or group of expres- 
sions. Implicit in both error criteria 
is the existence of requirement data r 
(e.g., loss equalizer requirements) to 
which computed values £ (e.g., loss of 
equalizer configuration) are to be matched 
by varying parameters x (e.g., element 
values of loss equalizer) . 



Figure 1. -Optimization Program Block 
Diagram 

Figure 1 outlines the steps required 
to solve typical engineering problems 
iteratively, and hence is also an outline 


for a computer program to automate the 
task. 

OPTIMIZATION METHODS 

The two general classes of optimiza- 
tion algorithms are the s lope -following 
and search techniques. The slope - 
following methods use a gradient or 
similar slope -determining function to 
find a new set of parameters which will 
reduce y. Search techniques are varied, 
but do not use gradient information. 

Slope -following techniques are often 
superior for functions with continuous 
derivatives, but are more time-consuming 
because of the large number of derivatives 
to be found. The search techniques, 
however, often work better with functions 
with discontinuous derivatives. In 
addition, constraints are often easier to 
include and, in most cases, a greater 
variety of optimization criteria can be 
incorporated. All optimization schemes 
except some random and grid search methods 
tend to reach the nearest local minimum 
of the function to be optimized. 

In the following sections the 
unrealistically simple function 

y = 16 x x 2 + (x 2 -4) 2 (4) 

will be used to illustrate various 
optimization techniques. 

SEARCH METHODS 
Grid and Random Search 

Grid search basically involves 
evaluating Equation (1) for many choices 
of x with each x^ taking on values 
throughout a prechosen range, as illus- 
trated on Figure 2. In random search the 
trial parameter values are chosen with 
some degree of randomness; in the illus- 
tration on Figure 3 the direction and 
distance of each move are chosen at 
random, and only those moves resulting 
in a lower y are retained. 

Pattern Searcn 

Pattern search techniques attempt to 
establish proper moves (succession of 
changes of x) by first making a series of 
exploratory maneuvers. After exploring 
the local terrain a larger move called a 
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SLOPE -FOLLOWING METHODS 


Steepest Descent Method 

The steepest descent method relies 
on the elementary calculus notion that 
the gradient of a function points in the 
direction of increasing functional values. 
Hence to minimize a function successive 
parameter sets can be calculated as 

*i + l = - k aCSi) • (5) 

Figure 8 illustrates a steepest 
descent path. The gradient partial 
derivatives also are used in many other 
algorithms, including the gradient partan, 
generalized Newton -Raphs on, and Fletcher - 
Powell methods. 

Gradient Partan 
13 14 

Gradient partan * is one of 
several versions of the method of parallel 
tangents. The method seeks to speed up 
the steepest descent approach by combining 
the results of gradient steps to choose a 
direction leading to a large improvement, 
as shown for the test problem on Figure 9. 
Considerable acceleration is often 
experienced. 

Generalized Newton -Raphson Techniques 

When x is close to the optimum most 
functions y are approximately quadratic. 

In the generalized Newton -Raphson tech- 
nique each X£ is replaced by xi+a£ Ax^. 

The partial derivative of the least 
squares error measure with respect to 
each a£ is then calculated and set equal 
to zero. The resulting equations are 
linear in a_, and can be solved. If the 
function g is nonlinear in x, as is 
usually the case, this procedure must be 
iterated. 

Fletcher -Powell Algorithm 

The Fletcher -Powell algorithm^ 
introduces a matrix H which is usually 
initially chosen to be a unit matrix. At 
each iteration H is varied, and it 
approaches the inverse of the matrix of — 
second derivatives. Instead of searching 
in the negative gradient direction for a 
minimum, a search is made along 

S = -(H) (Vy) . (6) 



A COMBINATION OPTIMIZATION PROGRAM 


The SUPROX Program 

The SUPROX (Successive apPROXimation) 
program^ combines the modified steepest 
descent algorithm with the generalized 
Newton -Raphson technique in a general 
program package.* The user must supply 
a routine to evaluate the function g, 
unless he can use one of the many routines 
in the program library file. Many speci- 
fied options are included, including 
limiting parameter variations to specified 
ranges, weighting requirements, and using 
requirement ranges. The program has been 
effective on many practical engineering 
problems, although no program has yet 
been found which will solve all problems. 


*The organization of the SUPROX program 
is due to Mr. C. L. Semmelman of Bell 
Telephone Laboratories, Holmdel, 

New Jersey. 
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RELATIVE DELAY REQUIRED (n sees) 


Sample Problem - Delay Equalizer Design 

A delay equalizer design problem is 
used to illustrate the use of SUPROX. An 
eight -section delay equalizer was required • 
as shown in Figure 10. The SUPROX func- 
tion evaluation subroutine works with 
ideal all-pass sections, and a subsequent 
program computes parasitic -corrected 
network element values from the idealized 
parameters . 







6 ADDITIONAL SECTIONS 




Figure 10 . -Eight -section Delay Equalizer 



-REQUIREMENTS 

AFTER 10 

/~ ITERATIONS 




TO 

FREQUENCY (MHz) 



FREQUENCY {MHz) 

Figure 11. -Delay Equalizer Requirements 
and Performance of Initial 
Design 


Without any attempt to design the 
equalizer, initial parameter values were 
chosen, producing the large deviation 
between required delay and network per- 
formance shown on Figure 11, 



FREQUENCY (MHz) 

Figures 12A,B,C. -Delay Errors at Various 
Stages of Optimization Process 



The error requirements and the orig- 
inal error are shown in Figure 12A, along 
with the results after four iterations 
when the modified steepest descent method 
has ceased producing significant improve- 
ments. The program then switched to the 
generalized Newton -Raphs on algorithm. 
Figures 12B and 12C illustrate the delay 
error after 10 and 47 iterations respec- 
tively. The requirements (+0.2 ns in the 
center region) were met with greater than 
a 2:1 margin remaining for physical real- 
ization errors. 

The performance of SUPROX can also 
be visualized by examining the match error 
y and the number of zero crossings of the 
delay error as a function of the number of 
iterations, as shown in Figure 13. 



Figure 13 . -Match Error and Number of Zero 
Crossings for Delay Example 


CONCLUSION 

The optimization techniques discussed 
in this paper have been applied to practi- 
cal engineering problems, but none is the 
desired panacea. The engineer should have 
access to programs with a variety of 
methods, and ideally should be able to 
interact with the computer during the 
optimization procedure. As new techniques 
and powerful time -sharing computer systems 
evolve the optimization procedure may 
again be dramatically altered. 
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Modular Structure 
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A progress report is presented on the first 
year of a cooperative effort to develop a modu- 
larly arranged circuit analysis program by 

(a) A group of seven universities with 
small to medium size computer faci- 
lities, each participant specializing in 
a module of the program. 

(b) Several industrial users of computer- 
programs with extensive experience, 
nearly all of whom are active in de- 
veloping their own in-house programs 
and who are willing first to trade know- 
how and secondly are willing to give 
serious thought to cooperative pro- 
gramming efforts. 

Present capabilities of NASAP are discus- 
sed including flowgraph construction and evalua- 
tion, input -output formulation and available 
subroutines. Strengths and weaknesses are pre- 
sented, together with plans for future develop- 
ment. 

NASAP DEVELOPMENT 

Introduction 

The Network Analysis for Systems Appli- 
cations Program (NASAP) has been developed by 
the NASA/Electronics’ Research Center as a two- 
fold project: 

(a) A research tool for the development of 
new concepts in automated circuit 
design. 

(b) An effective means of establishing 
qualification and standardization pro- 
cedures for future NASA computer 
programs. 

The first phase of the project has involved 
the definition and implementation of the symbolic 
oriented techniques, and an outline of their ef- 
fectiveness in a modularly arranged program. 
This work is well documented in the literature 
(1-12). This portion of NASAP will be discussed 
in detail and future work will be outlined. 


The NASAP program has been written in 
FORTRAN IV for the CDC 3600 computer. It has 
been designed in a modular arrangement, with 
the main part of the program designed to develop 
and evaluate a symbolic and numerical frequency- 
dependent signal flowgraph for a given network. 
The various modules perform specific applica- 
tions and extensions of the reduced flowgraph. 

This arrangement has provided an effective meth- 
od of dividing the program development among 
several universities and companies. 

The modular structure has also provided 
for expansion of the program to include additional 
subroutines or larger systems as the need arises. 
A more efficient use of computer time is achieved 
since the engineer selects which of several anal- 
yses or outputs he desires for a specific applica- 
tion, and bypasses those which he does not desire. 

The program has developed with seven ma- 
jor branches emanating from the flowgraph con- 
struction and evaluation program. These 
correspond to the areas of seven university grants 
in progress during FY67. Included are various 
display routines for the derived flowgraph, the 
addition of subroutines for non-linear elements 
Bode plot, transient response, sensitivity analy- 
sis, tolerance analysis and derivation and con- 
struction of an approximate flowgraph. In addi- 
tion, each university is responsible for qualifying 
the main portion of the program as well as its 
own subroutine on its own computer. This re- 
sults in a program which is qualified on a number 
of different computers. 

Maintenance Procedures 

To avoid costly duplication of effort for a 
large program (on failure diagnostics and exten- 
sion of scope) an atmosphere for a free exchange 
of program information is a prime requisite. 

For NASAP, program information and documen- 
tation is made available on the understanding that 
recipients will freely interchange programs, 
diagnostics, and additions. NASAP is a coopera- 
tive effort of about 20 users, whose development 
has fostered: 



(a) A modular structure of compatible 
subroutines, such that each subroutine 
can be improved, extended and elim- 
inted separately as needed. The main- 
tenance of each module is the assigned 
responsibility of one or more key 
participants with particular interests 
in these subroutines. 

(b) Qualification procedures, such as a 
sequence of sample problems with 
corresponding performance criteria; 
such as running times and round-off 
error for two or more machines. It . 
is also desirable to include a group of 
problems which reveal at what point 
the performance of the program be- 
comes marginal or unsatisfactory. 

Some preliminary results of this coopera- 
tive effort entailed in NASAP wiil now be dis- 
cussed. 

NASAP PROBLEM FORMULATION 

From a given network, NASAP constructs 
a corresponding block-diagram, often referred 
to as flow graph which relates voltages and cur- 
rents of the network elements. The block-dia- 
gram may have frequency dependent or time- 
dependent, linear and nonlinear, active and pas- 
sive elements. Active elements are modeled in 
the block-diagram either as voltage or as cur- 
rent sources which are controlled by the output 
of other sources. 

A passive element is modeled as a special 
case of an active element, in which source and 
control are of opposite type; namely, a voltage 
controlled current source (Y) or a current con- 
trolled voltage source (Z). Thus, for each pas- 
sive element a choice or dichotomy is necessary 
in modeling the network. This choice is con- 
strained by requiring that no voltage sources are 
in parallel and no current sources in series. 

The flowgraph approach, often referred to 
as the "dichotomous approach, " differs signifi- 
cantly from matrix oriented techniques. Like 
matrix techniques it serves to determine driving 
point and transfer functions, but it also provides 
valuable insight into problem formulation. An 
added feature is that network functions may be 
obtained in symbolic form. 

The NASAP program constructs a unique 
flowgraph for a given equivalent circuit. A net- 
work element in a flowgraph is represented as a 
transmittance G (S) relating two variables X and 
Y as in Figure 1. Evaluation of a transfer or 



PIG. 1 . Open System Representation, Y = GX 


driving point function is reduced to the determin- 
ation of the transmittance between an input and 
output node of the resultant flowgraph. For in- 
stance, for Figure 3 the problem is to determine 
the voltage gain. It is therefore necessary to 
determine the transmittance between the input 
node of the flowgraph (the voltage across element 
1) and the output node (the voltage across element 
9). In flowgraph notation this is as shown in 
Figure 1 where X = Vj and Y = Vo for Figure 3. 

For an actual circuit the corresponding flowgraph 
will be a many noded construction for which G(S) 
will not be at all obvious. The computer pro- 
gram evaluates this flowgraph by grouping al- 
gorithms which permit symbol manipulation within 
the computer. 

The flowgraph of Figure I is open, i.e. 
starting at the input node (X) there is no way of 
returning to that node in the direction of the trans- 
mittances of the flowgraph. Desired is a flow- 
graph which has only interdependent variables, 
such as Figure 2. To close the flowgraph a dummy 



FIG. 2. Closed System Representation, Y = GX, X = TY 

transmittance T is used. The program evaluates 
symbolically a closed system containing the un- 
known parameter 

T(S) = X (S)/Y (S) = 1/G(S). (1) 


It is instructive to interpret this result in terms 
of the circuit of Figure 3. To determine 
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FIG. 3. Example!: Two-Transistor Amplifier 

the voltage gain, the voltage of the input genera- 
tor, element 1, is made dependent upon the vol- 
tage across element 9, which is the output voltage. 
Hence, the desired transfer function V 0 /Vi can 
be calculated as 1/T(S). 

From the constructed block-diagram lin- 
ear, piece-wise, stochastic, or non-linear rou- 
tines can evaluate desired circuit criteria. 
Significant advances result from 

(a) Labyrinth and grouping algorithms, 
which circumvent combinatorial 
schemes used for tree enumeration in 
passive networks. 

(b) Tagging techniques, which utilize pro- 
perties of closed systems. For ex- 
ample, in symbolic calculations no 
distinction is made between known and 
unknown quantities. 

(c) Optimization procedures, which elim- 
inate from a circuit those design 
parameters which do not affect perfor- 
mance parameters more than a pre- 
assigned level of accuracy. Thus, the 
simplest possible model is displayed, 
given a set of design requirements. 

Network and transfer functions as well as 
sensitivity and generating functions are thus ob- 
tained symbolically and numerically in the fre- 
quency domain. State-space techniques yield the 
corresponding transient response in the time- 
domain or in the probability domain. 


Various display techniques can be used to 
construct Bode plot, pole-zero configuration and 
similar design graphs. 

\ 

NASAP INPUT -OUTPUT AND EVALUATION 
ROUTINES 

The tutorial description of NASAP to fol- 
low will outline the fundamental algorithms and 
describe in general terms the mechanics of the 
program, viz, : 

Coding of equivalent circuit (input 
statement) 

Construction and evaluation of flow- 
graph 

Function and operation of available 
modular subroutines 
Illustration by example 

The Equivalent Circuit 

The first step in preparing a circuit for 
computer analysis is the construction of a coded 
equivalent circuit from the given schematic. 

This serves one main purpose: the separation of 
all elements into two distinct groups (current 
generators coded "l", and voltage generators 
coded "0"). For passive elements, admittances 
and impedances are considered as voltage con- 
trolled current generators and current controlled 
voltage generators respectively. This dichotomy 
or separation of elements into two mutually ex- 
clusive groups is essential to the flowgraph con- 
struction by the computer. 

The elements of the coded schematic are 
numbered consecutively as are the independent 
voltage nodes, A group of elements selected as 
voltage generators is called the ’’tree”; those 
selected as current generators the "co-tree". 

A separation of the network elements into these 
two mutually exclusive groups, tree and co-tree, 
is made subject to certain constraints imposed 
by network physics: 

(1) Voltage sources in the circuit are en- 
coded as voltage generators. Thus in 
the coded equivalent circuit, voltage 
sources always form part of the tree. 
Conversely circuit current sources are 
encoded as current generators in the 
coded equivalent and hence never form 
part of the tree. 

(2) A passive element if contained in the 
tree is referred to as an impedance, 
(Z). If contained in the co-tree it is 
referred to as an admittance, (Y). 
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(3) When assigning elements to the tree 
and co-tree, voltage sources are never 
allowed in parallel or current sources 
in series. 

Coding of Equivalent Circuits for Computer 
Input 

The dichotomized equivalent circuit must 
be presented as input data to the program in the 
form of a matrix specifying circuit topology, 
frequency dependence, numerical values, and 
dichotomous structure. Table 1 gives this prob- 
lem statement for the circuit in Figure 3A. 

Each row corresponds to one circuit element 
numbered as in Figure 3B, and indicated in in- 
put E. 


Column A describes the vertex of origin of 
E, where the direction of the arrows in Figure 
3B specifies the assumed positive direction of 
current. Column B specifies the vertex termi- 
nation, Thus element 1 has entry A = 1, B = 2 
to specify the positive direction of current as 
flowing from node 1 to node 2. The inputs A and 
B thus describe the circuit topology. 

Column C indicates the dichotomy of vol- 
tage or current control to the program by the 
binary inputs ”0" for voltage and ”1” for current. 
Input D indicates the element which is the con- 
troller. This is the same as E for a passive 
element since a passive element performs its 
own control function between input and output. 
For the active controlled source E as 4, the cur- 
rent is controlled by element 3, hence D = 3. 
Input G serves to complete the dichotomous des- 
cription in the same binary manner as C. For 
example, from Table 1 and Figure 3B element 
3 is a passive voltage controlled current gen- 
erator, a resistor, hence is coded C = 0, D = 3, 
G = 1. 


As discussed in the second section (NASAP 
Problem Formulation), a closed flowgraph for 
the circuit is desired. This is specified in input 
H which identifies the unknown parameter, and 
indicates how the system is to be closed by the 
dummy transmittance T (S). For the example, 
the voltage gain is desired. Referring to (1), 
X(S) = V if Y(S) = V Q and T(S) = V^o- Hence 
1/T(S) represents the desired voltage gain. 
Element 1 is made dependent on element 9 for 
this calculation through input H(T) - 1, which 
represents 1/T(S) = G(S) for this problem. All 
other entries in column H are coded zero. 


The use of tagging parameters in the input 
serves to indicate quantities which require fre- 
quent and rapid identification. Column F lists 
the frequency dependence of the element E, tag- 
ging SF where F = -1, 0, +1 for electrical net- 
works. For example, a capacitor coded as an 
admittance has F = +1, since Y = SC . In this 
example all elements are frequency independent, 
hence F = 0 in ail cases. The program accepts 
higher powers of S, hence is useful for electro- 
mechanical and mechanical systems. 

Column K utilizes another tagging variable 
useful for sensitivity calculations. The entry 
K = 1 in the input for element 4 indicates to the 
program to calculate the sensitivity of the vol- 
tage gain to changes in the (B of transistor 1, 
since element 4 is the controlled current genera- 
tor of the transistor equivalent circuit. 

Care must be exercised in the specification 
of numerical inputs to scale the value of an ele- 
ment according to its dichotomous description. 
For example, a resistor coded as an admittance 
must have as its numerical input a value of con- 
ductance rather than resistance. Dependent gen- 
erators are coded with the value of the associated 
gain. 

Construction of Flowgraph 

Table 2 gives the output of the computer- 
constructed flowgraph in the form of three 
matrices V, W, and T, where N equals the num- 
ber of elements; here N = 9. These outputs rep- 
resent a judicious choice of 'the constraining 
equations for the closed system according to the 
tree chosen at the beginning of the problem form- 
ulation. 

The matrices V and W define the Kirchhoff 
voltage and current restraints for the network, 
hence represent equations for network inter- 
connections. The T matrix defines the set of 
differential equations for the network-intra- 
relationships for the elements. 

For example from Table 2, reading down 
column 1 of the W matrix: 

1(1) = -1(2) - 1(3) 

This is the Kirchhoff current constraint at node 

1 . 

It is to be emphasize^ that these outputs 
are only an intermediate step in the problem 
solution and hence serve only a tutorial function 
here. They are not visually necessary for the 
flowgraph construction and evaluation, since this 
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is performed within core, and these intermediate 
results are stored and used within core. 

Evaluation of Flow graph 

The basic concept in the evaluation of the 
fiowgraph is the first order loop L (l,n). A 
first order loop is defined as a connected se- 
quence of directed transmittances through the 
fiowgraph which touches no node more than once. 
All n — 1, N first-order loops are consecu- 
tively evaluated by labyrinth routines and are 
tabulated with their symbolic and numerical 
values. The value of a loop is the product of 
its constituent transmittances. 

A second order loop L(2, n) is defined as 
the product of two disjoint first order loops, i.e. 
those which have no common nodes. Similarly 
a loop of order m,L(m,n), consists of m dis- 
joint first order loops. These loops are succes- 
sively evaluated and presented in a similar 
manner to the first order loops. 

For a given closed system a constraint 
exists among the loops of the system fiowgraph. 

If there are N first order loops present, denote 
the sum of these loops by H(i) where 
N 

H(l) = 2 L(l,n) (2) 

n=l 

Let H(m) be the sum of all L(m,n). If the 
highest order loop present is of order M, then 
a M 

H = 2 (-l) m H(m) = 0 

m=0 (3) 

where H(0) = 1 and (3) is defined to be the 
topology equation constraint for closed systems 
which equals zero. H is thus one plus the signed 
sum of all loops present in the fiowgraph. 

The topology equation takes the form of a 
linear function of the fictitious closing para- 
meter T, (t') 

H = H(T) + T H = 0 (4) 

where H(T')is that part of the equation devoid 
of T 

and H(T*) is that part deprived of T 

(since T appears only linearly, it can 
easily be factored out by the computer 
using the tagging technique described) 

Equation (4) can then be solved for the desired 
function G = 1/T yielding 

G = 1/T = -H(T')/H(T) 


which is the desired frequency-dependent transfer 
function. 

Table 3 gives the computer fiowgraph eval- 
uation for example 1. The columns L(i) and H(i) 
give the number and order of the loops starting 
with H(0) numbered 0; E(i) defines which of the 
nine transmittances are contained in the loops. 

The columns T, K, and S indicate the tags for 
each loop which are formed as the algebraic sum 
of the tags of the constituent elements. Finally, 
the numerical value of each loop is listed, with 
H(0) = 1. The transfer function is then formed 
according to the procedure outlined above. The 
computer does this by forming separate signed 
sums according to (3) of those loops which are 
tagged with T = 0 or T = 1 since H(T) » H(T = 0) 
and H(T’) = H(T = 1). The quotient is then formed 
according to (5) to yield Vg = 9.00. Table 4 gives 
the calculation of the sensitivity of the voltage 
gain with respect to 0, which is discussed in more 
detail just below. For now, just note that this 
sensitivity = 1, showing that 0 has a direct effect 
on the voltage gain of the circuit. 

DESIGN APPLICATIONS 
Sensitivity Analysis 

The fiowgraph approach is particularly 
useful in the calculation of sensitivity functions. 
We define the sensitivity of a performance 
criterion P with respect to a parameter Q as 

P d log P „ dP/P 
S = d log Q dQ/Q 
Q 

The program then calculates the symbolic and 
numerical sensitivity function defined above. 

The details of the algorithm are not im- 
portant here, but it is important to note that the 
technique eliminates the need for numerical 
schemes to take the required derivatives. These 
are accomplished instead by grouping algorithms 
for the symbolic topology equation. (13) 

For the problem of example 1, the sensi- 
tivity of the voltage gain to variations in 0 is 
required. It can be shown that 

Vg ZHfPj) 2H(K=0) 

S 0 x = ZH(vgj = SH(T=0) 

which can be obtained from Table 3. 

H(0 1) and H(Vg) are the parts of the 


(5) 
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topology equation devoid of 0 j and Vg. A com- 
plete description of the sensitivity analysis por- 
tion of NASAP is available. 


Transient and Frequency Response 

A documented subroutine (14) is available 
to obtain the transient response of a network 
from the calculated transfer function. This sub- 
routine produces the response of the network to 
impulse, step, ramp and sinusoidal inputs . An 
arbitrary input may be specified as an array of 
numbers. The outputs consist of a graphic dis- 
play of the response as well as a table of num- 
bers showing the response of the system at 
equally spaced intervals of time. 

A plotting subroutine (Bode Plot) Is also 
available to display the frequency dependent 
transfer function over a specified range of fre- 
quency giving both magnitude and phase. 

Example 2: L-C Filter 


To further illustrate the techniques used, 
and as an example of a circuit with frequency 
dependent elements, consider the L-C filter of 
Figure 4A. The dichotomized schematic is 
shown in Figure 4B while the flowgraph drawn 
from the computer output in Figure 4C is shown 
only as a visual aid. It is desired to find the 
input admittance of this circuit, and its sensi- 
tivity to variations in L2. 




FIG. 4. Example 2: L-C Filter 


The problem statement or input to the pro- 
gram is shown in Table 5, which has been formed 
according to the rules developed in the subsec- 
tion on Coding of Equivalent Circuits for Com- 
puter Input. Note particularly row 1 of this 
array for which C = 0, D = 2, G = 1, and H = 1. 
This indicates that the unknown transmittance is 
one which relates the voltage generated by ele- 
ment 2 and the current generated by element 1: 
the input admittance. Column F of this matrix 
contains the appropriate entries for the frequency 
dependent elements. For example, element 3 
is a capacitor coded as a voltage controlled cur- 
rent generator (C=0, G=l) or admittance, has 
frequency dependence of Si hence F (E=3) = +1. 
The K = 1 in row 2 indicates that the sensitivity 
analysis is to be performed with respect to L2. 

Table 6 gives the flowgraph evaluation as 
before. The transfer function must now be form- 
ed according to the tags T=0 or 1, but also group- 
ing those terms having the same S tag to form 
the appropriate coefficients of terms with like 
powers of S. An additional output here is the 
transfer function normalized as to the power of 
S and the coefficient of the highest power of S in 
both numerator and denominator polynomials. 
Table 7 gives the sensitivity analysis, also is a 
quotient of polynomials according to the previous 
subsection on Transient and Frequency Response. 

It should be noted that these are not the 
only outputs available. Additional subroutines 
factor all polynomials of interest, form a 
transient response and a Bode plot. Planned 
sub-routines, as well as further work on the 
above routines are discussed in the concluding 
section. 


CONCLUSIONS 

The problem statement of NASAP seems 
at first appearance somewhat cumbersome to 
use, although with a little practice it becomes 
quite easy. 

The choice of the tree at the very begin- 
ning of the coding seems to be quite critical to 
the running time of the problem. There is no 
theoretical information available on the best 
(least computing time) choice of a tree, but ex- 
perience has shown that a closely connected tree 
(many voltage generators emanating from one 
node) has a shorter running time than one which 
is loosely connected. For most practical net- 
works, the ground node naturally becomes this 
node. 
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The program at present handles a linear, 
time- invariant circuit with a maximum of 40 
elements. Most practical networks can be 
adequately approximated by a linear equivalent 
circuit only in fairly narrow ranges of operation. 
For this reason it is planned to extend the pro- 
gram to piecewise linear and non-linear networks. 

Flow graph techniques may be used effec- 
tively to reduce the complexity of circuit models. 
An optimization technique is to be implemented 
as a subprogram which finds and eliminates all 
components which contribute to a desired result 
less than a pre -assigned value. 


The NASAP program has received wide- 
spread attention because of its symbolically- 
oriented algorithms and its inherently logical 
output capabilities. In this respect it has ful- 


filled the initial reasons for its development as 
a research tool. 

Grants and contracts are now in force to 
extend and develop NASAP into a useful design 
tool. The capacity is being extended to 120 
elements. The input formats are being rewritten 
to provide automatic generation of the input 
matrix. But, since the construction of the matrix 
puts the designer in close touch with the circuit 
he is analyzing, it has been decided to leave this 
automatic generation as an option. 

A complete plotting package is being 
written for display of the various outputs at the 
user's request. A tolerance analysis and a 
sensitivity matrix giving the sensitivity of trans- 
fer functions of interest to all components now 
available as subroutines, are being added mod- 
ular ly to the main program. 
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TABLE 1 


Problem Statement for Two -Transistor Amplifier 


A 

B 

C 

D 

E 

F 

G 

H 

K 

Numerical 

1 

2 

0 

9 

1 

0 

0 

1 

0 

.100E+1 

1 

2 

0 

2 

2 

0 

1 

0 

0 

. 100E-2 

1 

2 

0 

3 

3 

0 

i 

0 

0 

. 100E-1 

3 

2 

1 

3 

4 

0 

1 

0 

1 

. 100E+2 

3 

2 

1 

5 

5 

0 

0 

0 

0 

. 100E+4 

3 

2 

0 

6 

6 

0 

1 

0 

0 

.100E-1 

4 

2 

1 

6 

7 

0 

1 

0 

0 

. lOOE-f-2 

4 

2 

d 

8 

8 

0 

1 

0 

0 

. 100E-2 

4 

2 

i 

9 

9 

0 

0 

0 

0 

. 100E+2 


TABLE 2 


Coiiiputer-Output-Flowgraph Construction 



TABLE 3 


Fiowgraph Evaluation for Two-Stage Amplifier 


L(i) 

H(i) 

m 

TKS + 

Numerical 

0 

0 


000 + 

.100E+1 

1 

1 

1 34567 9 

110 - 

.100E+5 

2 

1 

56 

000 + 

. 100E+2 

3 

1 

89 

000 + 

. 100E-1 

4 

2 

56 89 

000 + 

.lOOE+O 


Transfer Function = .900E+1 
Normalized Transfer Function = .900E+1 


TABLE 4 


Calculation of Sensitivity Function for 
Two-Stage Amplifier 


m 

Numerical 

S 

0 

.100E+1 

0 

H(Pi) 2 

.100E+2 

0 

3 

. 100E-1 

0 

4 

. 100E+0 

0 


Sum = . U1E+2 


0 

.100E+1 

0 

H(Vg) 2 

.100E+2 

0 

3 

. 100E-1 

0 

4 

.100E+0 . 

0 


Sum = . 111E+2 


Sensitivity = .100E+1 



TABLE 5 


Problem Statement For L-C Filter 


A 

B 

C 

D 

E 

F 

G 

H 

K 

Numerical 

1 

2 

0 

2 

1 

0 

1 

i 

0 

. 100E+1 

1 

2 

1 

2 

2 

1 

0 

0 

i 

.666E+1 

1 

2 

0 

3 

3 

1 

1 

0 

0 

.120E+1 

1 

3 

0 

4 

4 

-1 

i 

0 

0 

. 200E-1 

3 

2 

1 

5 

5 

-1 

0 

0 

0 

. 120E+3 
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TABLE 6 2. 

Flowgraph Evaluation For 
L-C Filter 


3. 


4. 


5. 


6 . 


TABLE 7 

7. 

Sensitivity of Voltage Gain of L-C Filter 
to Variations in L 2 


8 . 


9. 


10 . 


1. R. M. Carpenter and W. W. Happ, 

Algorithms for the Dichotomous Repre- 11. 

sentation of Macrocircuits, Proceedings . 

Fourth Annual Allerton Conference on 
Circuit and System Theory, October, 1966. 
Accepted for publication by IEEE Trans - 
actions on Systems. Science, and Cyber- 
netics. for Aug. , 1967. 



L(i) 

Numerical S 


0 

.100E+1 0 

H(L 2 ) 

4 

. 240E+1 -2 

Sum = 

, 100E+1 + (,240E+1)S -2 


0 

. 100E+1 0 

H(Zin) 

2 

. 880E+0 2 


3 

. 133E-1 0 


4 

.240E+1 -2 


6 

. 192E+1 0 

Sum = (,880E+0)S 2 
Sensitivity = 1.14— 

+ 2.9 + (,240E+1)S" 2 
S 2 + 2.40 


s' 1 

+ 3.31S 2 +2.72 
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SUMMARY 


N 67-22634 


The basic features, of ten automated nonlinear 
transient circuit analysis programs are listed and 
the values of each of the features discussed from 
a potential user’s viewpoint. The features c can- 
pa red for each program include the types of analy- 
sis performed, the types of elements handled, 
built-in models, the computers on which the pro- 
grams are operational, documentation, program 
availability, user convenience, special features, 
and mathematical solution formulation. 


I. Capability of the Program: This section 

was concerned with the types of analysis per- 
formed (ac, dc, transient), kinds and number of 
elements handled, and built-in model capability. 

II. Computer Compatibility: This section 

asked on which computers the program was opera-* 
tional, in what language the program was written, 
the size memory required, and questions concerning 
the status and documentation of the program. 


INTRODUCTION 

This paper compares the major features of 
several digital computer programs for nonlinear 
transient circuit analysis. The features of the 
programs arc presented from the user's viewpoint 
and include discussions of the capabilities, limi- 
tations, and availability of each program. 

There are in existence a number of transient 
circuit analysis computer programs. Most of these 
programs can be applied to Transient Radiation 
Effects on Electronics (TREE) problems. Since 
there is no current, comparative documentation on 
the capability, limitations, and availability of 
these programs, the Defense Atomic Support Agency 
(DASA) requested the Air Force Weapons Laboratory 
(AFWL) review the programs in existence. 

APPROACH 

This survey complements earlier surveys by 
Wirth 1 in 1964 which briefly compared PREDICT, 
NET-1, CIRCUS, MISSAP, and ECAP; by Dickhaut 2 in 
1965 which compared PREDICT, NET-1, and CIRCUS; 
and by Pritchard^ in 1965 which compared TRAC with 
PREDICT, NET-1, and CIRCUS. Programs reviewed 
previously have been modified, new programs have 
been developed, and several other industry- and 
university-developed programs have become of gen- 
eral interest through recent listings in tech- 
nical journals. ^>5 

Information for this survey was primarily 
ga thered from a questionnaire which was sent to 
sixteen different government laboratories, cor- 
porations, and universities which were known to 
have developed sophisticated circuit analysis 
programs . 

The questionnaire was formulated to gather 
information in a concise form about all aspects of 
the program of interest to a potential user of 
the program. The questionnaire consisted of seven 
sections : 


III, Radiation Effects: This section was con- 
cerned with special features of a program for 
handling radiation effects problems. 

IV. Use, Convenience, and Flexibility Fac- 
tors: This section was concerned with input and 

output formats, restrictions on circuit topology 
specification, and general features such as re- 
run options. 

V. Mathematical Details: This section was 

concerned with the types of solution formulation 
and integration routines used in the program. 

VI. Availability: This section was concerned 

with the availability of the program and the pro- 
cedure for obtaining the program. 

VII. Additional Information: This section 

asked for a user's manual and a sample problem 
run for each program. 

RESULTS 

From the information obtained from the ques- 
tionnaires,* ten programs are categorized as auto- 
mated nonlinear transient circuit analysis pro- 
grams. These programs are automated in that they 
formulate the circuit equations from topological 
and component data and perfoim the solution cal- 
culations to provide a transient or time history 
analysis of the circuit. The programs are also 
powerful enough to include nonlinear elements 
(capacitors, resistors, inductors, or voltage 
or current sources) used to model active elec- 
tronic devices. 

The names of the ten programs included in 
this survey are given in Table I along with the 
names of the originating organization and 


*See Acknowledgments section of this paper and 
references 6 through 13 . 
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government sponsor, A comparison of the major 
features and capabilities of the programs is shown 
in Tables II through IV. 

DISCUSSION 

The program features are compared in Tables 
II through IV. The following is a user-oriented 
discussion of these features. As seen in Table II 
all programs except PREDICT provide a steady-state 
dc solution. Each program uses an iterative tech- 
nique, usually, a modified Newton -Raphs on to obtain 
the dc solution. Mils provides an efficient means 
of obtaining the initial conditions required for a 
subsequent transient analysis. Two of the pro- 
grams, the General Network Analysis program and 
the MISSAP I also provide nominal ac solutions and 
frequency or Fourier analysis. 

One of the first concerns of a transient 
analysis is the initial conditions of the circuit. 
For all programs, except PREDICT, the initial con- 
ditions for the transient analysis can be automat- 
ically calculated by the dc solution portion of 
the program. For PREDICT the dc solution can be 
Obtained by a separate "power supply-turn-on" 
transient run. All of the programs except CIRCUS 
and the General Network Analysis program provide 
for entering user -supplied initial conditions. 

This feature is particularly useful for circuits 
which the dc iterative solution fails to converge 
and the circuit initial condition must be approxi- 
mated from measurements or separate transient runs 

An important feature of each program is the 
type of elements or components which may be 
entered. All of the programs readily accept con- 
stant valued resistors, capacitors, inductors, and 
constant and time-varying voltage sources. All 
except NET-1 also include constant and time- 
varying current sources. All except TAG and TRAC 
provide for entering inductive coupling through 
mutual inductance. These elements along with 
built-in models for active devices provide an 
adequate tool for describing most standard dis- 
cret component transistorized circuits. For 
applications such as modeling new semiconductor 
devices, more flexibility is required. For these 
applications it is important to be able to define 
functionally -variable elements and enter these in- 
to the program. Several of the programs provide 
for this by allowing variable elements defined by 
a table (F 2I , equation (F 3 ), or subroutine (F J, 
as shown in Table II. 

An indication of the maximum size of the cir- 
cuit which may be entered is shown in the next two 
colums of Table II. This varies from 30 nodes for 
the General Network Analysis program and 60 ele- 
ments for the Oklahoma State Analysis program, to 
300 nodes for the SCEPTRE program and 600 elements 
for the NET-1R program. For some of the programs, 
additional limits are placed on the number of each 
type of element, while others limit only the total 
number of elements or nodes. 

The CIRCUS, General Network Analysis, NET-1R, 
and TRAC programs include fixed built-in models 


for active devices. These models are considered 
fixed models in that they cannot readily be 
changed by a user. The built-in transistor and 
diode models of these programs are nonlinear tran- 
sient Ebers-Moll models and are essentially the 
same for each of these programs. 

MISSAP III built-in models are defined in 
subroutines which can be changed by a user. Other 
programs, PREDICT, SCEPTRE, STRAP, and TAG do not 
include built-in models. For these programs ac- 
tive devices are entered by user -defined equiva- 
lent circuit models either as part of an overall 
circuit, or in the case of SCEPTRE, from a user- 
defined model library tape. 

The last column on Table II indicates whether 
models or model parameter data may be qtored and 
called from a model library tape. The SCEPTRE- 
stored model feature is unique in that the stored 
model is user defined and may have up to 25 exter- 
nal terminals. 

Table III indicates the computers on which 
each of the programs are operational, the language 
it is written in, the extent of the documentation, 
and the availability of each program. 

If a program is capable of solving a particu- 
lar problem and is operational on an available 
computer, whether or not the program is used may be 
determined by how difficult the program i3 to uBe. 
All of the programs have been designed to be used 
by engineers without requiring the services of a 
programmer. All of the programs use an engineer- 
oriented input language. In general, the less 
restrictive, easier to use programs provide for a 
free-field input as indicated in Table IV. 

An example free-field SCEPTRE input is shown 
in Figure 1. The input for this circuit would 
look similar for CIRCUS, with commas replacing 
the dash and equal signs, and for NET-1 with 
blanks replacing the commas, dashes, and equal 
signs. The input for PREDICT would be similar 
except the transistor would have to be specified 
ns equivalent circuit elements and defining equa- 
tions. STRAP would also appear similar except 
variable element values would be entered In a 
subroutine. 

An example of the free-field subroutine for- 
mat of MISSAP III is shown in Figure 2. 

TAG also uses a free-field subroutine format 
which is flexible, but less convenient to use. 

TRAC uses a fixed field format for element speci- 
fication plus the use of auxiliary FORTRAN equa- 
tions for defining nonstandard elements or outputs. 

An economically important feature is a save 
and continue feature that permits the user to 
save the results of a transient run in a form 
that permits the run to be continued from that 
point. This is particularly valuable for analyz- 
ing large circuits which may require several 
minutes of computer time for a few microseconds 
of problem time and where the circuit recovery 
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time cannot be estimated accurately. All programs 
except General Network Analysis\and MISSAP III 
have this feature. © 


R1=1KJ1 
JWV 



CIRCUIT DESCRIPTION 

ELEMENTS 

Rl, 1 -2=1 

R2, 2-4= 19 

*3, 5 - 3 = 1.5 

KL, 1-5-10 

ET, 4 - l = 10 

Tl, 2 - 1 - 3 = MODEL 2N914A 

OUTPUTS 

VR3, TO, VCXT1, PLOT 
RUN CONTROLS 
STOPTIME - 500 
RUN INITIAL CONDITIONS 
END 


Figure 1, Example of SCEPTRE Input Format 


Another convenient feature of most of the 
programs is one that permits the circuit analysis 
to be rerun automatically, with circuit parameter 
changes (element value, transient forcing func- 
tion, etc.) by requiring only specification of the 
changes . 

For all the programs, transient output is a 
time history of the circuit response. Some of the 
programs such as CIRCUS, General Network Analysis, 
and NET-1R essentially provide only node voltages 
and semiconductor junction currents and voltages 
ns output. PREDICT provides only element currents 
and voltages as output. Other -programs such as 
SCEPTRE, TAG, and TRAC also provide for additional 
user -defined outputs which are combinations of 
circuit variables. Using this feature, quantities 
such as element power dissipation and voltage 
between arbitrary points in a circuit are readily 
available for output. MISSAP III provides a 
unique method for specifying output by inserting 
voltmeters and ammeters in the circuit and then 
printing these "meter readings" as output. 

Most of the programs provide for plotted re- 
sults as well as printed listings as shown in the 
"plots" column of Table IV. Some of the programs 
such as CIRCUS and SCEPTRE provide plot informa- 
tion which can be processed by user installation 
supplied plot routines to obtain additional plot- 
ted results . 



CALL MISSAP ( 6 l, 40, 1000) 

CALL E (2, 6 , 3, 1.5, 0, 1.0E-3, 1000., 0) 

CALL V (2, 6 , 1, 1, 1.0E-3) 

CALL T ( 2 , 4, 2 , 15 1.0E-6, 0.98, 1 . 0 E- 7 , 39 . 5 ) 
CALL R (4, 6 , 330. O) 

CALL C (4, 6 , 0.05E-6. 0.7, 0) 

CALL R (15, 8 , I 8 OO.O) 

CALL DC ( 8 , 6 , 12.0) 

CALL V (15, 8 , 2, 1, 1.0E-3) 

CALL QRUN (0, 5-OE-3, 5-OE-4) 

RETURN 

END 

Figure 2, Example of MISSAP III Input Format 

The next four columns of Table IV indicate the 
solution termination options available for each 
program. Most of the programs provide for solu- 
tion termination when the problem response time 
limit specification is reached, when a specified 
computer machine time limit is reached, and when 
the solution time step increments decrease below 
a minimum limit. MISSAP III, SCEPTRE, and TAG 
also provide for termination when circuit vari- 
ables exceed user-specified limits. 

One of the major uses of these programs is 
providing circuit nuclear radiation response 
calculations. All but three of the programs, 

MISSAP III, Oklahoma State Systems Analysis, and 
TAG were either designed or have been adapted 
specifically for radiation effects applications. 

Some of the programs provide directly for 
radiation effects analysis as shown in Thble IV 
by including photocurrent generators in the built- 
in semiconductor device models. These programs 
also usually provide a simple format for entering 
radiation effects data. These programs also pro- 
vide for entering or calculating radiation in- 
duced device model parameter changes. 

The remaining columns of Tfeble IV indicate 
the basic approach used for the equation formula- 
tion and solution by each of the programs. For 
the transient solution either a state variable 
(capacitor voltage and inductor current) or nodal 
(admittance matrix) formulation is used. Gener- 
ally, an explicit numerical integration routine 
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is used with the state variable f omuls t ion. Gen- 
erally, for the nodal f emulation, implicit inte- 
gration by including recursive difference approxi- 
mations in the admittance matrix is used for 
solution with an accompanying iteration solution 
technique for nonlinearities and dependent sources 
at each solution time Increment. 

There are other important aspects of these 
programs that have not been included in this paper. 
One is the relative solution efficiency of the 
various programs. Previous comparisons have show 
differences of up to a factor of ten in the comput- 
er time required for different programs to cal- 
culate the same circuit response on a similar 
computer. Although efficiency is important, its 
accurate evaluation requires that all programs be 
run on the same computer for a family of circuit 
problems . 

CONCLUSIONS 

The programs reviewed were quite similar in 
capability of analyzing most transistorized cir- 
cuits. The entering of most circuit data is con- 
venient and quite similar for CIRCUS, NET-1R, and 
SCEPTRE, with TAG and TRAC probably the least 
convenient. 

Considerable differences in the flexibility 
of the programs were noted. CIRCUS, General Net- 
work Analysis, and NET-1R programs are fairly 
inflexible in that only fixed value elements are 
permitted outside of the built-in models. The 
other programs are capable of processing func- 
tionally variable elements and so are useful for 
performing modeling work. However, the ease of 
entering variable elements varies considerably be- 
tween programs. Only PREDICT, SCEPTRE, and STRAP 
make provisions for describing variable elements in 
simple engineer -oriented language. 
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SABLE X— AUTOMATED NONLINEAR TRANSIENT CIRCUIT 
ANALYSIS PROGRAMS, ORIGINATORS, AND SPONSORS 


Program 

Originator 

Sponsor 

CIRCUS 

(Circuit Simulator) 

Boeing Company 


General Network Analysis 
Program 

Lockheed 
Sunnyvale, Calif 


MISSAP III 

(Michigan State Systems 
Analysis Program; 

Michigan State University 
East Lansing, Michigan 


NET-1R (Network Analysis 
Program. Radiation 
Version; 

Los Alamos Scientific 
Laboratory and Braddock, 
Dunn, & McDonald, Inc 

LASL (AEC) 
HDL, US Army 

Oklahoma State Systems 
Analysis Program 

Oklahoma State University 
Stillwater, Oklahoma 


PREDICT (Prediction of 
Radiation Effects by 
Digital Computer 

IBM Corporation 
Owe go, New York 

AFWL, USAF 

SCEPTRE (System for Cir- 
cuit Evaluation and Pre- 
diction of Transient 
Radiation Effects) 

IBM Corporation 
Owego, New York 

AFWL, USAF 

STRAP (Simplified Tran- 
sient Radiation Analysis 
Program) 

Douglas Aircraft Co. 
Santa Monica, Calif 


TAG (Transient Analysis 
Generator) 

Jet Propulsion Laboratory 
Pasadena, Calif 

NASA 

TRAC (Transient Radiation 
Ana lys i s by Computer ) 

Autonetics 
Anaheim, Calif 
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TABLE II: MAJOR FEATURES AND CAPABILITIES 
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* MISSAP I Element Values 

r* in state vector C = Constant value. F2 = Function of? circuit variables defined by a table 

;* Ideal transformer F = User-defined variable value F3 = Function of circuit variables defined by an equation 

Fi •= Function of time Ft, = Function of circuit variables defined by a subroutine 


















TABLE m— COMPUTER COMPATEBILTM, DOCUMENTATION , AND AVAILABILITY 
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TABLE IV: CONVENIENCE, SPECIAL FEATURES, AND MABMATICAL DETAILS 
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1620 CARDS - BCD CODE 

The circuit cards that have a BCD Code can 
still be used with System 360 if a card with the 
following characters is inserted in front of the 
circuit deck. 

■♦‘•BCD Columns 1-i*. 

All the following circuit decks will be as- 
sumed in the BCD mode unless a new mode card of 
the following form is encountered. 

-♦‘-EBCDIC Columns 1-7 

This card returns the system to the 
EBCDIC 360 mode. 


DC - AC TRANSIENT OPTIONS 

The following calculation can be made in 
all three types of analysis* 


CODE 


1. 

NODE VOLTAGES 

NV 

or Voltages 

2. 

ELEMENT CURRENTS 

CA 

or Currents 

3. 

ELEMENT VOLTAGES 

CV 


h. 

BRANCH CURRENTS 

BA 


5. 

BRANCH VOLTAGES 

BV 


6. 

ELEMENT POWER LOSS 

BP 



1 ERROR 


N 67-22635 


2 i B 

' i E is a voltage source in branch j 

TeT J 


* 1 L is a current source in branch j 
k* ^ ®i GM is transconductance in branch j 

3'gm, 3 


The sensitivity coefficients are defined 
as the change in node voltage for a one-percent 
change in the branch parameter. The four sensi- 
tivity coefficients available in DC analysis are 
the following* 


i 

- node 

* E i 



100 

}Ei 



100 




100 

JEi 



100 


The control card 1 ERROR - P is available in 
all three analyses* 

A summation is made of the currents at all 
nodes and if the SUM of the unbalances is greater 
than the 1 ERROR value the unbalances and the 
branch currents are printed. 


DC ANALYSIS OPTIONS 


II* Wors t Case Solu tion 

Tolerance data must be supplied in the 
following manner for a worst case solution* 

R = 1000. (900., 1100.) or R = 1000.(.l) 

E = 20. (19., a.) or E = 20. (.05) 

The worst case maximum solution is de- 
fined as the sum at the k^ 11 node as 


I* Partial Derivatives and Sensitivity Co- 
efficients 

The partial derivatives refer to the rate 
of change of voltage at a particular node 
with respect to a circuit parameter in a 
particular branch. The four partial de- 
rivatives available in DC analysis are the 
following* 


i - node j - branch 


1. <) E i r^ i s a Resistor in branch j 


hi A * 


Cjc - nominal vali 


where all the terms A P L are positive 
JPi 


A P^ - Tolerance of the parameter 

If a worst case minimum solution is de- 
sired, all the terms are negative. 


After the nominal solution has been cal- 
culated, parameter values are chosen for 
the worst case minimum solution. For the 
minimum solution, the minimum value of a 
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Miscellaneous 


parameter is chosen if its partial derivative is 
positive. If the partial derivative is negative, 
the maximum value is used. A new solution for 
node voltage is made using the new parameter 
values. The partial derivative is calculated 
again, and if a sign change from the nominal case 
is encountered, a warning message is printed out. 
This warning message indicates that a true worst 
case minimum has not been found. 

Next, the worst case maximum solution is 
obtained. The procedure is the same as that of 
the minimum solution except that if the partial 
derivative is positive, the maximum parameter 
value is substituted and if the partial de- 
rivative is negative, the minimum parameter 
value is substituted. Again, if a sign change 
in the partial derivative is encountered, a 
warning message is printed out and the true worst 
case maximum has not been found. 

III. Standard Deviation 

The calculation of standard deviations 
requires that the user supply the i 3 standard 
deviation values of the circuit parameters. 
These values are entered the same way that 
tolerance data is entered. An assumption is 
made that the circuit output variables are 
linearly related to the parameters, so that 
the partial derivatives are nearly constant 
over the range between maximum and minimum 
values of the parameters. 


DC-AC OPTIONS 
Parametric Variation 

Parameters can be varied only in DC and AC 
analysis. Far example, to change a branch that 
contained a resistor 

B1 N(0,l), R=10 

another two cards of the following form would be 
placed after the last branch card. 

B 13 N(5,0), R= 500 LAST BRANCH CARD. 

MODIFY 

B 1 R = 5(3) 20 

PRINT, VOLTAGES 
END 

The resistor in branch 1 would be incre- 
mented from 5 to 20 ohms in 3 steps (5 -ft per 
step). The node voltages would be printed out 
for each increment. 

A single value could have been substituted 
in the following manner 
MODIFY 

B1 R = 20 


The MISCELLANEOUS output block indicator 
can be used in DC and AC analysis only. 

When MISCELLANEOUS is encountered in a DC 
analysis print command, the nodal admittance 
matrix, equivalent current vector, and nodal im- 
pedance matrix will be printed. 

In an AC analysis, only the nodal admit- 
tance matrix and the equivalent current vector 
will be printed if MISCELLANEOUS is encountered 
in a print coirmand* 


AC OPTIONS 
Frequency Variation 

Frequency variation can be either linear or 
geometric. For linear frequency variation, the 
following two cards are entered 

MODIFY 

FREQ = PI (+P2)P3 

PI - starting frequency 

P3 - end frequency 

P2 - number of steps 

For geometric frequency variation, the 
following two cards are entered 

MODIFY 

FREQs PI (P2)P3 

PI and P3 are the same as above 
?2 - multiplier 


TRANSIENT ANALYSIS OPTIONS 

Time Dependent Voltage and Current Sources 

The card containing the source information 
goes behind the branch that the source is lo- 
cated in. The card format is the following 

Xnn (K), PO, FI, P2, ... Pn. 

X can be E or I nn = branch number 

K= TIKE STEFS the source value is specified. 
PO, PI, P2, ...Pn = source values at selected 
time intervals.. 
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Equilibrium 


Bsriodic So urce 

Xnn P(K), PO, PI, P2, ... Pn 
\ 

indicates a periodic source 

Sinusordal So urce 

Xnn SIN (TP), VI, VO, TO 



If an EQUILIBRIUM control card is placed 
after the print statement, only the steady state 
solution will be printed out. 


2 ERROR 

The 2 ERRORS p card determines the time 
within which switch actuation takes place. The 
actuation will take place in P x TIME SIEP of 
the TIME STEP in which actuation occurs. 


3 ERROR 

The card 3 ERROR - P can reduce the time 
step after an initial condition solution. If 
P>1, no reduction of the time step takes place. 
If P<1, then the first time step immediately 
after the initial condition solution, is sub- 
divided into four parts. 


CIRCUIT I 



CIRCUIT 2 
LOW PASS FILTER 



CIRCUIT 3 



SWITCH IS OPENEO 
AT T = 0 
V,2 (0-) =2.5 


-1.0.9- 




CIRCUIT 4 

TRANSISTOR SWEEP CIRCUIT 




- 110 - 



imB3707E 










112 ' 





/ 

























115 . 











■118 



• 11 9 - 





120 
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the application of computers in the engineering design process. 
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PROGRAM DESIGN FOR SCIENTIFIC COMPUTING 


by 

K. E. WITTMER 

Bell Telephone Laboratories, Incorporated 
North Andover, Massachusetts 


ABSTRACT 


The design of computer programs is a very 
essential part of computing and often it 
does not receive enough attention to guar- 
antee economical computer usage. 

Scientific computing is reviewed from the 
standpoint of the computer user. 

Program design considerations are presented 
with emphasis on the development of appli- 
cation programs. 

Important program design principles are 
demonstrated by means of a user-oriented, 
special-purpose network analysis program. 


Because of review complications. Dr. Wittmer 
was unable to submit his complete paper in 
time for publication. 
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